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Preface 


The purpose of this study was to investigate the nonlinear behavior of a 
simple po'-ared 1 .'.sting hypersonic vehicle flying in a near circular orbit 
above a spherical nonrotating Earth with gradients in atmospheric density 
and pressure and an inverse square law for gravity. The vehicle is 
constrained to fly in a vertical plane so only longitudinal motion is 
modeled. Bifurcation analysis, utilizing the AUTO software package, was 
used to conduct this study. A simple five-state model with three 
different thrust laws was used to describe an unaugmented vehicle whose 
geometric and aerodynamic characteristics follow those of the literature. 
A parameter represented a body flap deflection (5^) was used to conduct 
one set of bifurcation sweeps for each thrust law. Then a second set of 
bifurcation sweeps for each thrust law was obtained using a parameter 
representing a throttle (6T) which scaled the value of the thrust. 
Secondary parameters representing simple feedback gains, were subsequently 
added. 

I wish to extend my sincerest thanks to my thesis advisor, Capt Jim 
Planeaux, for his patient and caring nature. His guidance and insight 
were invaluable. I would also like to thank the members of my thesis 
committee, Dr. Brad Liebst and Major Curtis Mracek for their comments 
while reviewing this document. Finally, but most importantly, I would like 
to t..^nk my dearest friend and love, Cynthia, for always being there for 
me and shouldering the responsibilities of our family during my time at 
AFIT. 

Eric E. Fox 


ii 







Table of Contents 


Page 

Preface .....ii 

List of Figures .....v 

List of Tables ...viii 

Notation .....ix 

Abstract .xii 

I. Introduction .1 

Introductory Discussion .1 

Summary of Previous Studies . 2 

Outline of Analysis .9 

Equations of Motion .12 

Aerodynamic Forces and Moment Coefficients ..16 

Thrust Laws . 18 

Vehicle Characteristics .20 

II. Introduction to Bifurcation Analysis and Stability .22 

Equilibrium Solutions and Equilibrium Points .22 

Simple Nonlinear Behavior .23 

Limit Cycles (Orbits) .28 

Two Parameter Continuation .31 

III. Bifurcation Analysis of Longitudinal Dynamics.33 

Body Flap Parameter (5uj) Variation .. 36 

Constant Thrust Rocket Case . 39 

Air-Breathing Engine Case .53 

Throttle Parameter (5T) Variation .61 

Constant Thrust Rocket Case .61 

Air-Breathing Engine Case .63 

IV. Model Stabilization .69 

Simple Feedback Options .69 

V. Conclusions & Recommendations .72 

Appendix A: FORTRAN Listing of the Required AuTO Subroutines .76 

Appendix B: Standard Atmospheric Approximations for Density 

and Pressure .64 


iii 




































Page 

Bibliography .87 

Vita ...89 










List of Figures 


Figure Page 

1. Frequencies near Resonance Altitude (15:19) .6 

2. Axis System Relative to Inertial Reference .11 

3. Comparison of C^/Cp with Etkin's Work .21 

4. Pendulum for Example Problem .25 

5. Bifurcation Diagram for Example Problem .27 

6. Limit of Periodic Motion about a Solution Branch .29 

7. Limit Cycle or Orbit in Phase Space (12:63) .29 

8. Curve.- of Hopf Points for Two Parameter Continuation 

with Pitch Rate Gain and Body Flap Deflection .32 

9. Bifurcation Diagrams for Body Flap Sweep from 

100 kft : Constant Thrust Rocket .40 

10. Bifurcation Diagrams for Body Flap Sweep from 

300 kft : Constant Thrust Rocket ..40 

11. Left Branch Limit Cycles for 5^ = -56.6° and 
-41.8° Body Flap Sweep from 300 kft : Constant 

Thrust Rocket .43 

12. Right Branch Limit Cycles for 5^ = 56.73° and 
42.25° Body Flap Sweep from 300 Kft : Constant 

Thrust Rocket . 43 

13. Right Branch Limit Cycles in Phase Plane for 

= 56.73° and 42.25° Sweep from 300 kft 
Constant Thrust Rocket .44 

14. Left Branch Limit Cycles in Phase Plane for 

= -56.64° and -41.8° Sweep from 300 kft 
Constant Thrust Rocket .44 

15. Bifurcation Diagrams for Body Flap Sweep from 

400 kft : Constant Thrust Rocket .46 

16. a Bifurcation Diagram for Body Flap Sweep from 

400 kft : Constant Thrust Rocket...47 

17. h Bifurcation Diagram for Body Flap Sweep from 

400 kft : Constant Thrust Rocket .48 


v 





















Figure 


Page 


18. Limit Cycles for Periodic Branch 1 (5 bf = -6.826° & -2.83°) 


Body Flap Sweep from 400 kft : Constant Thrust Rocket .50 

19. Limit Cycles in Phase Plane, Periodic Branch 1 

(5faf = -6.826° & -2.89°) : Const Thrust Rocket ...50 

20. Limit Cycles, Periodic Branch 2, (5 bf = 8.738° & 3.57°) 

Body Flap Sweep from 400 kft : Const Thrust Rocket .51 

21. Limit Cycle in Phase Plane, Periodic Branch 2, 

(5 bf = 8.738° & 3.57°), Sweep from 400 kft : Const Thrust 
Rocxet . 51 

22. Growth of Nonlinear Waveform, Periodic Branch 2, 

a Bifurcation Diagram, Body Flap Sweep from 400 kft .52 

23. Bifurcation Diagram for Body Flap Sweep from 100 kft 

Air-Breathing Engine .54 

24. Partial Body Flap Sweep from 100 kft for the 

Air-Breathing Engine Case .55 

25. Limit Cycles for 5 b£ = -52.77° & 8 b£ = -20.54°, 

Body Flap Sweep from 100 kft : Air-Breathing Engine .56 

26. Limit Cycles in Phase Plane (5 b{ = -52.77° & -23.54°) 

Body Flap Sweep from 100 kft : Air-Breathing Engine .57 

27. Collection of Bifurcation Diagrams for Body Flap Sweep 

from 400 kft : Air-Breathing Engine .58 

28 Limit Cycles over one Period 5 bf = -0.0007° and 0.0118° 

Body Flap Sweep from 400 kft : Air-Breathing Engine Case .59 

29. Limit Cycles in Phase Plane for 5v £ = -0.0007° and 
0.0118° Body Flap Sweep from 400 Kft : Air-Breathing 

Engine . 59 

30. Limit Cycles, Left Branch, (5 b£ = -0.0046° & -0.011°) 

Body Flap Sweep from 400 kft : Air-Breathing Case ..60 

31. Limit Cycles in Phase Plage, Left Branch 

(5 b j =-0.0046° & -0.011°), Sweep from 400 kft : 

Air-Breathing Case .60 

32. Collection of Bifurcation Diagrams for Throttle Sweep 

from 100 kft : Constant Thrust Rocket . 62 


vi 


















Figure Page 

33. Collection of Bifurcation Diagrams for Throttle Sweep 

from 400 kft : Constant Thrust Rocket .62 

34. Collection of Bifurcation Diagrams for Throttle Sweep 

from 100 kft for Air-Breathing Engine Case .64 

35. Expanded View of Region near Hopf Point for 5T Sweep 

from 100 kft for Air Breathing Engines .65 

36. Maximum Limit Cycles for the 5T Sweep from 100 kft for 

Air-Breathing Engine Case .65 

37. Limit Cycles over One Period for 6T = 19.06 & 19.07 

Throttle Sweep from 100 kft : Air-Breathing Engine .66 

38. Expanded View of a Limit Cycle, 5T = 19.063 

Throttle Sweep from 100 kft : Air-Breathing Engine .67 

39. Limit Cycles in Phase Plane for 6T = 19.063 

Air-Breathing Engine Case from 100 kft .67 

40. Limit Cycles in Phase Plane for 5T = 19.078 

Air-Breathing Engine Case from 100 kft .68 

41. Comparison of Circular Orbital Period of Limit 

Cycles for Throttle Sweep from 100 kft : Air-Breathing .68 

42. Pitch Rate Feedback Loop for Model Stabilization .70 

43. Movement of Hopf Bifurcation given Pitch Rate 

Feedback to the body flap. Air-Breather form 100 kft .71 

44. Comparison of Calculated Density to Standard 

Atmosphere ...86 

45. Comparison of Calculated Pressure with Standard 

Atmosphere .86 






















List of Tables 


T 

* 



Table • Page 

1. Summary of some Feedback Techniques 

used by Stengel (13:470)...8 

2. Coefficients for Density Polynomial .8 


viii 





Notation 


a : angle of attack (radians) 

8 bf : body flap deflection (degrees) 

5 b £ C : commanded body flap deflection (degrees) 

5T : throttle parameter 

A : perturbation quantity 

y : flight path angle (radians) 

X : parameter for example problem (=F/mg) 

p : longitude (radians), also eigenvalue 

6 r pitch angle (radians) 

p : density (slug/cu ft) 

: angular velocity of body frame to inertial frame (radians/sec) 

: angular velocity of Earth (radians/sec) 

u 71 : angular velocity of vehicle carried frame to 

inertial frame (radians/sec) 

tt v: : angular velocity of wind frame to inertial frame (radians/sec) 

An : area of exhaust nozzle (ft ) 

a : acceleration (ft - sec *) 

C D : drag coefficient 

Cpg : basic drag coefficient 

Cp a : partial derivative of drag w.r.t. angle of attack (1/radian) 

Cp : lift coefficient 

Cpo : lift coefficient at zero angle of attack 

Cpg : partial derivative of lift w.r.t. angle of attack (1/radian) 

C, : pitching moment coefficient 

E 

Cjjj : basic pitching moment 
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: partial derivative of pitching moment w.r.t. 

body flap deflection (1/deg) 

: partial derivative of pitching moment w.r.t. 

angle of attack (1/radian) 

: partial derivative of pitching moment w.r.t. 

pitch rate (sec/radian) 

: drag (lb) 

: force generated by constant torque (F=T/r) in example problem (lb) 
: gravity as a function of radius from Earth's center (ft - sec’^) 

: altitude (ft) 

: (I, - I,)/I, 

: radius of gyration in pitch (ft) 

: pitch rate feedback gain 
: characteristic length (vehicle length) (ft) 

: lift (lb) 

: mass (slugs) 

: sum of moments about vehicles center of mass (ft - lb) 

: Mach number 
: ambient pressure (psf) 

: nozzle exit pressure (psf) 

: pitch rate of vehicle relative to Earth (radian/sec) 

: pitch rate relative to the Earth at starting altitude (radian/sec) 
: geocentric radius (ft), also length of rod in example problem (ft) 
: area of lifting surface (ft*) 

: thrust (lb), Period (seconds), torque (ft-lb) 

: velocity (ft/sec) 

: rocket nozzle exhaust velocity (ft/sec) 

: weight (lb) 


x 







X r : partial derivative of the longitudinal force w.r.t. radius 

X u : partial derivative of the longitudinal forces w.r.t. velocity 

in the & x direction 

() e : equilibriuni value; also value for example problem (i.e. © e ) 

{) 0 ; values at initial starting equilibrium solution 

(°) : derivative with respect to time 

(I x ,Ij,I x ) : moments of inertial about body axes (slug - ft 15 ) 

(£, j, £) : unit vectors for the inertial frame 

(£ x ,B y ,£ e ) : unit vectors for the body fixed frame 

($ x , ~Qy, ^ z ) : unit vectors for the vehicle carried frame 

( tf x , $ y , ti z ) unit vectors for the wind axis frame 
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Abstract 

Bifurcation analysis was used to investigate the nonlinear behavior 
of a simple powered lifting hypersonic vehicle in circular orbit about a 
spherical nonrotating Earth with gradients in atmospheric density and 
pressure and an inverse square law for gravity. Vehicle motion is 
constrained to a vertical plane so only longitudinal dynamics were 
modeled. Bifurcation analysis was conducted using the AUTO software 
package. A simple five-state model with three different thrust laws was 
derived to describe an unaugraented vehicle whose geometric and aerodynamic 
characteristics follow those of the literature. A parameter representing 
a body flap deflection {5 b j) was used to conduct one set of bifurcation 
sweeps for each thrust law. A second set of bifurcation sweeps for each 
thrust law was obtained using a parameter representing a throttle (5T) 
which scaled the thrust. Secondary parameters representing simple 
feedback gains were subsequently added. Results were surprising for a 
simple system with basically linear aerodynamics. Periodic branches 
arising from the loss of pitch stability or associated with a "resonance 
altitude" are routinely found with significant amplitude, and periods on 
the order of an elliptical orbit's period for a given geocentric radius. 
Rotational states generally had sub-oscillations of greater frequency. 
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BIFURCATION ANALYSIS OF THE LONGITUDINAL DYNAMICS 
OF A SIMPLE POWERED LIFTING HYPERSONIC VEHICLE 


I. Introduction 


Introductory Discussion 

In the past interest in hypersonic vehicle dynamics has concentrated 
on assuring stable reentry and return to a few specific points on the 
earth. The need for maneuverability was limited. The renewed interest in 
lifting hypersonic vehicle dynamics and design; brought about by the 
Trans-Atmospheric Vehicle projects, Boost Glide Vehicles, the National 
Aerospace Plane and other hypersonic lifting vehicles designed for 
improved maneuverability and greater versatility; has generated a need to 
better understand and anticipate their possible nonlinear dynamic effects. 
Specifically, the ability to predict nonlinear dynamic responses, periodic 
equilibrium states and other dynamic phenomena for a representative 
hypersonic vehicle is growing in importance. Previous work by Etkin (4), 
Berry (2), Vihn (15) and others demonstrate some interesting behavior of 
the longitudinal dynamics for hypersonic vehicles due to the variation of 
atmospheric density, gravity and Mach number with altitude. Bifurcation 
and continuation analyses have been used successfully to examine the 
nonlinear behavior of fighter aircraft in a variety of configurations. It 
was felt this type of analysis would yield insight into the nonlinear 
behavior of a hypersonic vehicle as well. 
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The purpose of this thesis was to explore the nonlinear dynamic 
responses of a hypersonic vehicle and using a more global technique to 
investigate these effects. In addition, it is hoped the application of 
bifurcation analysis techniques to the highly nonlinear hypersonic regime 
would help extend the basic techniques available for further analysis of 
hypersonic vehicles. 

Summary of Previous Studies 

Several papers have been presented over the last four decades that 
impact directly on the study of the longitudinal dynamics of a hypersonic 
vehicle flying in an atmosphere that contains gradients in density and the 
effects of curvature of the flight path. Host of these previous works 
built in some way upon the work presented in 1950 by Neumark (9), which 
was then extended to a lifting vehicle in orbital flight by Etkin (4) in 
1961. The results of these later works have served to enhance the 
material originally found in these two landmark papers for lifting 
vehicles. Having said this one should note that some correction of 
Etkin's observations regarding the behavior of the phugoid and pitching 
mode characteristics are found in the work by Vihn and Dobrzelecki (15) 
and verified in a later paper by Harkopoulos, et al (7) as well as this 
author's most recent work. 

In his paper Neumark details the motivation for bis work which was 
based on several of the very first studies of the behavior of airplanes in 
steep angled dives. He is one of the original writers on the subject of 
the effect of density gradients on airplanes having published his first 
work on the subject in 1931; the earliest being in 1929. His paper 
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published in 1950 was the first published work in which the longitudinal 
equations of motion were cast in the form of a quintic? having added an 
equation to describe the change in altitude with time. This form of the 
longitudinal equations gives rise to a fifth, real root (eigenvalue) which 
gives an indication of the vehicles ability to hold a fixed equilibrium 
altitude. The results of his study demonstrate the increasing affect the 
density gradient has on the longitudinal dynamics as the speed of the 
vehicle approaches Mach one. Neumark found that the principal effect of 
the density gradient is on the phugoid mode? he states that the short 
period (pitching) mode is unaffected. Increasing density, increases the 
phugoid frequency thus shortening the period. The effect on phugoid 
damping was not clear having been complicated by compressibility effects 
at speeds above H *1.4 . He concluded the height mode would have a very 
long time constant and may be either a subsidence or divergence and has 
importance only for hypersonic flight or flight at constant altitude for 
long periods of time (9:325). 

Etkin's classic of 1961 extends Neumark's work to the truly 
hypersonic case and includes, necessarily, the mathematical modifications 
to account for the curvature of the undisturbed flight path and the 
variation of gravity with altitude. In his analysis the longitudinal 
equations of motion for flight in a vertical plane about a nonrotating 
Earth whose atmosphere is at rest are linearized and the behavior of the 
vehicle subject to small disturbances about an equilibrium is examined. 
In addition, he presents results from numerical solutions to the nonlinear 
equations and does a comparison with the linear approximation. Of note in 
the equations of motion is the addition of the torque about the vehicle's 
center of mass due to the small variation of gravity acting on a body at 
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very high altitudes (above 500,000 ft) where the pitch damping is 
negligible. In this realm the gravity torque generates the dominant 
moment acting on the vehicle and for a standard vehicle configuration 
whose longitudinal axis is nearly aligned with the flight path this effect 
is destabilizing. 

Etkin examines four basic cases using the same (steady reference) 

lift coefficient (C t = 0.05 [rad'*]). 

Case A Constant thrust rocket, full set of equations 
Case B Air-breathing engine (T « p), full set of equations 

Case C Approximate equations (i.e. no density gradient) 

Case D Constant thrust rocket with qj = 0 

where qg is the equilibrium value of the pitch rate relative to the Earth. 

Etkin determined that the effects of varying density and gravity 

with altitude and the effects due to the Earth's curvature and the thrust 

law have significant impact on the phugoid mode and the stability of the 

height mode, but insignificant effect on the pitching (short period) mode 

except at very high altitude where the pitch damping is overcome by the 

gravity torque. In addition, Etkin found that above 400,000 ft the period 

of the pitching and phugoid modes approached each other and he asserts 

that they become equal, after which the phugoid tends toward the orbital 

period and the pitching mode tend toward infinity. In this altitude range 

he demonstrated a dynamic coupling between the two modes and when nearly 

equal all relation to two classical modes breaks down with substantial 

pitching motion in the phugoid mode. Finally when the two modes are 

exactly equal he determines the system to be unstable. For the height 

mode Etkin found that it represents "a spiral, proceeding away from the 

reference orbit." He also noted an interesting variation in the way the 

speed changes with altitude above and below 350,000 ft. Above this 
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altitude as the altitude is decreased the speed increases whereas below 
this altitude the opposite is the case (4:787-738). 

In the work by Vihn and Dobrzelecki (15) an "analytic study of the 
longitudinal dynamics of a thrusting, lifting orbital vehicle in a nearly 
circular orbit" is conducted. The basic set of five equations used to 
describe the longitudinal motion of a vehicle in orbit about a spherical 
Earth were used. A strictly linear analysis as well as analysis including 
second order terms in the Taylor series expansion of the atmospheric 
density were used to develop explicit relationships to describe the 
orbital motion. Also developed were analytic expressions for the period 
and damping of the "angle of attack" (pitching) mode. As with Etkin they 
observed an altitude where the velocity-altitude relationship inverts. 
They went on to develop an expression based on vehicle characteristics 
that defines the altitude where this "inversion" takes place. Finally 
they found the trend at high altitude of the linearized phugoid or long- 
period mode and angle of attack (pitching) mode tend to become nearly 
equal in frequency, period and damping, then diverge. (Figure 1) Similar 
behavior for the very same equations was found earlier by Etkin (4:785- 
788) where he concluded the two modes "crossover" and the phugoid period 
tended to the orbital period while the period of the pitching (short 
period) mode tended toward infinity. At the point of "crossing" Etkin 
concluded the dynamic system would be unstable (4:787). 

Stengel also found the same basic trends in the three longitudinal 
modes however his work looked more closely at the stability questions and 
dealt at some length with various techniques to provide altitude stability 
for a vehicle in supersonic cruising flight (13). In his work he uses the 
linearized equations for longitudinal motion, characterized by the 
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Figure 1. Frequencies near Resonance 
Altitude (15:19) 

perturbation variables Au (forward velocity), Aa (angle of attack), A© 
(pitch attitude), and Ah (altitude), to study the interrelations of 
motions these variable characterize. He then used this information to 
test how various combinations of feedback and control would affect 
altitude stability. In addition to developing analytical transfer 
functions he conducted numerical studies and summarizes the effectiveness 
of the various techniques proposed. Some of his results were tabulated in 
his work and are reproduced on the following page in Table 1. 

In a later work by Berry (2), he examined the effect on the "long- 
period dynamics" of a vehicle of similar characteristics to previous 
authors but included an "advanced air-breather" with a more complicated 
thrust law which is more representative of a true hypersonic vehicle like 
the National Aerospace Plane. He concluded that the height mode stability 
and long-period damping were strongly affected by the slope of thrust with 
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Mach number (6:255-257). While all the trends indicated in his work are 
valid they are not in the strictest sense complete. It is the variation 
of the difference between thrust and drag with height and Mach number that 
determines stability. These more complete relations are identified in 
works by Markopoulos, et al (7:285) and Stengel (13:468,472). Berry 
further examined the effectiveness of various simple feedback schemes 
involving only the pitch control surface to stabilize the long-period 
and/or height modes. Some results for feeding back combinations of pitch 
attitude, forward velocity and altitude to the pitch control surface for 
the rocket and "advanced air-breather” were presented (2:256-257). 

Most recent is the work by Markopoulos, Mease and Vihn (7) where the 
linearized equations of motion are used to examire the thrust law effects 
on the longitudinal dynamics of an aerospace vehicle flying at hypersonic 
speeds. Their work demonstrates the dependence of the height mode 
stability and phugoid damping on the way the longitudinal force varies 
with both altitude and speed. In addition they confirm the results of the 
previous study by Vihn and Dobrzelecki (15) for the high altitude trend of 
the phugoid or long-period mode and angle of attack (pitching) mode. 

Of special interest in the work by Markopoulos, et al, is the 
characterization of expected height mode stability and phugoid damping 
over the plane of all thrust possibilities. The plane is defined by two 
parameters, specifically the variation in the longitudinal force with 
height (X r ) and velocity (X u ). The correlation of the points relating to 
earlier work by Etkin (4) are in excellent agreement. They went on to 
conclude that it is actually "the partial derivatives of the difference 
between thrust and drag with respect to speed and altitude that plays the 
key role in determining the stability of the translational dynamics (7:287)." 


7 




TABLE 1 


SOME FEEDBACK EFFECTS FOR STENGEL'S STANDARD CASE WITH AUGMENTED SHORT 
PERIOD (FEEDBACK IS NEGATIVE UNLESS DENOTED BY (+)) 


Feedback Variable 


and Control 

Height Mode 

Phugoid Mode 

Attitude to: 

thrust 

SS 

SI 

lift (-) 

I 

S 

lift (+) 

SS 

I 

moment (-) 

S 

I 

moment (+) 

SI 

S 

Pitch angle to: 

thrust 

N 

I 

lift (-) 

I 

SS 

lift (+) 

N 

SI 

moment (-) 

I 

ss a 

moment (+) 

N 

SI 

Forward Velocity to: 

thrust 

SS 

ss a 

lift (-) 

SI 

SS 

lift (+) 

N 

SI 

moment 

SS 

ss a 

Angle of attack to: 

thrust 

N 

SS 

lift (-) 

SI 

SI a 

lift (+) 

N 

SS 

moment (-) 

? 

SS b 

moment (+) 

? 


S = Stability SS = Strong Stability 

I = Instability SI = Strong instability 
N = Neutral stability 
3 With Limited Travel. 

b Conditional Stability Stengel (13:470) 

mmmmammmmmmmmmmmmmmmmmmmmmmmmmmmmmaMummmmmmmmmmamamBmmmmHmmmmmmmmmmmmam 


This same relationship was discussed by Stengel in his article on 
"Altitude Stability for Supersonic Cruising Flight" (13:468,472). 
Markopoulos, et al, further concluded after numerical simulation of their 
full and reduced order mathematical models that over the plane of all 
engine possibilities 

"if thrust increases faster than drag vith respect to speed at 
least one of the translational modes (height or phugoid) vill 
be unstable. Increasing the partial derivative of the differ- 
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ej ice between thrust and drag with respect to altitude has a 
destabilizing effect on the height mode and a stabilizing 
effect on the phugoid (7:287)." 

Finally, they concluded that to first order, the period of the phugoid 
mode as well as all characteristics of the pitching mode are independent 
of the thrust law (7:287). 

Outline of Analysis 

In this thesis the dynamic behavior of a powered lifting hypersonic 
vehicle in nearly circular orbital flight about the center of mass of a 
spherical nonrotating planet (specifically the Earth) whose atmosphere 
contains gradients in density and pressure and whose gravitational field 
follows the inverse square law is examined. 

Throughout this thesis emphasis is given to the leading order 
aerodynamic behavior and simplifying assumptions to this end are brought 
to the readers attention as required. In order to focus the scope of this 
work it is assumed the flow is inviscid therefore the effects of high 
temperature gas flows are neglected. This assumption is consistent with 
general longitudinal analysis found routinely in the literature and allows 
use of simple Newtonian impact theory as the basis for the aerodynamics. 

To begin this study the reader should have a good mental image of 
the problem being analyzed, and a well developed understanding of the 
equations used to describe the translation and rotation of a body flying 
a great circle about a spherical planet. To facilitate this the basic 
equations for a vehicle flying in an atmosphere at rest relative to a 
nonrotating spherical planet (tt ? =0) are derived. The first step in 
analyzing this problem is to identify an inertial reference frame. Then, 
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three advantageous frames of reference relative to the inertial frame are 
introduced from which a set of equations describing the forces and moments 
acting on the body of interest are developed. It is common in trajectory 
analysis to use a wind axis system as shown in Figure 2 where the positive 
x-axis is parallel to, but opposite, the relative wind. In this way the 
aerodynamic forces are cleanly defined and the velocity vector has a 
single non-zero component. For the analysis of angular momentum the body- 
fixed axes are used, also shown in Figure 2 as b x and b z , thus the moments 
of inertia are time invariant and for a fixed mass and mass distribution, 
as is the case here, the moments of. inertia are constant. The vehicle 
axes indicated in Figure 2 by (v)_ , is used as a convenient intermediate 
frame between the body or wind axes and the inertial axes. 

As with linear analysis, the bifurcation analysis begins at a known 
equilibrium point, but rather than linearizing the equations and looking 
at small disturbances about this point, a continuation method is used to 
solve for the flow of equilibrium solutions (specifically the pseudo-arc 
length technique resident to AUTO; the software package used for 
continuation and bifurcation analysis in this study) (3:12-16; 12:116). 
From the path of equilibrium solutions (or stationary points) bifurcating 
solutions, or in other words, additional solution paths are located and 
explored. Of the various types of possible solution branches special 
interest is given to branches obtained subsequent to a nonhyperbolic or 
degenerate point (12:18-20). These often give rise to branches of 
periodic solutions where motion, such as periodic oscillations develop. 
On these branches the dynamic behavior comes to life. Many of these 
concepts are clarified in section II where the nature of bifurcation 
analysis is discussed. 


10 









11 









Equations of Motion 


The author is indebted to Planeaux (10), HcRuer, et al (8:204,220) 
and Etkin (5:104,148) for background and guidance for the following 
development. The reader is encouraged to review the two texts for details 
on the development of the following system of equations. To begin, the 
basic assumptions on which the equations of motion are based must be 
stated: 

1. The Earth is an inertial reference frame. 

2. The vehicle is a rigid body. 

3. The vehicle mass and mass distribution are constant. 

4. The vehicle is symmetric about the x-z plane. 

5. The body fixed axes are aligned with the principal axis of the 
vehicle. 

In addition to the simplifications resulting from the assumptions 
above, the terms associated with motion in the horizontal plane are 
neglected leaving the system as shown above in Figure 2. 

Looking first at the angular and kinematic relations, by inspection the 
following angular rate of the three reference frames relative to the 
inertial frame are given as: 

= gj = (-£ + $) j » (6 + -£) £ y 


= (-£) J = "|i tfy 


( 2 ) 


12 





( 3 ) 


s± wi = (-£ + Y) J = (Y - A) ^y 


where: (i, j,£) 

(£ xl S yl £ s ) 

(ti x ,fi yl & 2 ) 

Y 

6 


i* vi 

q 


unit vectors of inertial frame 

unit vectors of body frame 

unit vectors of wind frame 

time rate change of longitude (rad/sec) 

time rate change of flight path angle 

(rad/sec) 

time rate change of pitch angle (rad/sec) 
angular velocity of body to inertial frame 
(rad/sec) 

angular velocity of vehicle to inertial 
frame (rad/sec) 

angular velocity of wind to inertial frame 
(rad/sec) 

pitch rate of the vehicle relative to the 
Earth (rad/sec) 


The radius (r) is the distance from the center of mass of the Earth to the 
vehicles center of mass and written as a vector in the vehicle frame is 
given as: 


r = 




(4) 


where: = unit vectors of the vehicle frame 

r = geocentric radius (ft) 

From eqn(4) the velocity can be written by differentiating in the vehicle 

frame: 


V = £ = -£V Z + iii vi xr 
= -i $ 2 + 


where: y 

V 


flight path angle (rad) 
velocity (ft/sec) 


13 



However the velocity can also be shown to be: 

V= V(cos (y) V x - sin(y) Q g ) (6) 

By combining eqn(5) and eqn(6) the following scalar kinematic relations 
are obtained: 

i = Vsin(y) (7) 

|i =-pCos(y) (8) 

where: y = flight path angle (rad) 

r = radius (ft) 

V = velocity (ft/sec) 

From eqn(l) the equation for the time rate change of pitch angle can be 
obtained by solving for 6 as follows: 

6 = g + |i (9) 

Following directly from the assumptions 3 through 5 the time rate change 
of the pitch velocity of the vehicle is given as (4:145): 

$ = -y- (10) 

where: I y = moment of inertia about the y-axis of the body (slug-ft J ) 

= sum of moments about the vehicle^ center of mass (ft-lb) 

Equations (7,9 and 10) make up three of the five equations of 

motion. The remaining two equations fall out of the force balance which 

is dealt with next. 
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Summing forces acting at the vehicles center of mass yields: 


F - m& = mg(z) V z + tS x - L$ z - D$ x 

= [-mg(r) sin(y) + Tcos(a) - D]# x 

+ Dn£rcos(y) - rsin(o) - L] 


(11) 


where: a = acceleration <ft-sec" 2 ) 

D = drag (lb) 

g(r) = gravity as a function of geocentric radius (ft-sec -2 ) 

L = lift (lb) 

m = mass (slug) 

r = geocentric radius (ft) 

T = thrust (lb) 

Now velocity in the wind frame is written as V=Vti x . The acceleration 
in the wind frame is given by: 


V = a = V# x - (y-(x) Vti z 


( 12 ) 


Since F=ma eqn(12) and eqn(ll) can be set equal, after multiplying 
eqn(12) by the mass m, and upon separating into scalar components yields 
the following two equations. 


mV = Tcos(a) - D - mgsin(y) 


(13) 


-mtf - |i) V - -Tsin(a) - L + mg cos(y) (14) 


Utilizing the relation 0 = y + a and eqn(9) the following expression is 
obtained: 


g-d = 


(15) 
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Substituting eqn(15) into the left hand side of eqn(14) and solving ford 
yields the final equation for the set of five equations of motion. 

& = g - -4sin(«) - -4 + J^ilicostY) (16) 

mV mV r 

Equations (7,9,10,13 and 16) comprise the set of five dynamic and 
kinematic equations for this analysis. Together with the following 
expressions for the aerodynamic forces and moments, and the thrust 
equations they comprise the complete set of equations required to conduct 
this study. The five equations of motion are reprinted below for 
convenience: 


mV = Tcos(a) - D - mgs in(y) 

(13) 

& = G - -^sin(a) - — + “4-cos(y) 
mV mV i 

(16) 

6 = g + A 

(9) 

• L y 

(10) 

£ = Vsin(y) 

(7) 


Aerodynamic Forces and Moment Coefficients . As with linear analysis 
the standard forms cf the forces and the moment due to aerodynamics will 
be used, and are listed below. Notice however, the term on the right hand 
side of eqn(19). This term is the moment about the center of mass of a 
satellite in a gravitational field and as found by Etkin is a significant 
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factor at high altitudes (4:783; 11:21) 



D = 



(17) 


L = 

-- f P V“SC t 


(18) 

*-■§ 

P v 2 sic b 

- A £L (j - j ) 
2 r x 3 

sin (20) 

(19) 


ii 

Cp 0 + 9D6 


(20) 


C L * 

CLO + CLa a 


(21) 

^JD = CjO 0 

+ 

* (g - ji) + 

C Mr ibf 

(22) 


where: 


p = atmospheric density (slug/ft 3 ) 

C D = nondimensional drag coefficient 

C D0 = basic drag coefficient 

C Da = partial derivative of C p w.r.t. alpha (1/rad) 

= nondimensional lift coefficient 
C^q = basic lift coefficient 

C La = partial derivative of C t w.r.t. alpha (1/rad) 

Cj = nondimensional aerodynamic moment coefficient 

Cj, = partial derivative of C | w.r.t. pitch rate (sec/rad) 

Cjg = partial derivative of C a w.r.t. alpha (1/rad) 

C,t bf = partial derivative of C. w.r.t. body flap deflection 

(1/deg) 

g(r) = gravitational acceleration (ft/sec 2 ) 

1 = characteristic length (vehicle length) (ft) 

S = area of lifting surface (ft 2 ) 


The final term in eqn(22) represents the contribution of the body flap 
deflection to the moment coefficient and is used as a standard pitch 


control surface. The density is calculated using one of two analytic 
expressions depending on the altitude. The specifics of how the density 
is calculated as well as a brief discussion of the development of the 


analytic expressions is found in Appendix 2. 


17 







Thrust Laws . Three basic ideal thrust laws are used in this study. 
The equations representing these thrust laws are detailed below. 

constant thrust rocket 

T = -|p 0 vfcCi* (23) 

where: p c = atmospheric density at starting altitude (slug/ft 3 ) 

Vq = velocity at starting altitude (ft/sec) 

Cj,q = drag coefficient at starting altitude 

Notice for the constant thrust case the thrust is fixed at the values of 

the drag for the starting altitude (i.e. the altitude where the 

bifurcation sweeps starts). This is required since to start the 

continuation method an equilibrium solution must be provided as a first 

step. In all cases the equilibrium solution has a = 9 = 0 radians. This 

requires the thrust to equal the drag for equilibrium. 

variable thrust rocket ( 6 : 356 ) 

T = itiV^ + <P 0 - P a ) A d (24) 

As stated above the starting equilibrium point with a = 0 = 0 radians 
requires the thrust to equal drag when the continuation method is begun. 
This requires that the mass flow be determined by setting the thrust equal 
to the drag at the starting altitude, therefore the mass flow rate of the 
exhaust is fixed at the following equation. Note it is assumed the mass 
flow is sufficiently small relative to the mass of the vehicle as to be 
negligible. This assumption is fairly good at high altitude but is very 
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poor at altitudes below about 200,000 ft. 


ih = 


-f Po - (P, - P.) A d 


(25) 


where: 



exhaust nozzle area (ft^) 

mass flow rate of exhaust (slug/sec) 

pressure at exhaust nozzle exit plane (psf) 

ambient air pressure (psf) 

velocity of exhausted mass (ft/sec) 


ideal turbojet (9»33 2 ) 


T = 



(26) 


As with the constant thrust rocket, thrust for the ideal turbojet must 
equal drag at the starting altitude since the continuation method requires 
an initial "starting" equilibrium solution and a = 8 = 0 was taken as 
the values of these states at the equilibrium point. Therefore the value 
of Tq is set by the following relation: 

T 0 = -|Po VoSCoo (27) 

Note the exponent "x" can be used to change the way the thrust varies with 
altitude. For the standard turbojet, x equals one (x=l). 

It should be noted that in all cases the thrust can be varied with 
in the subroutine CONST through a Thrust Scaling Parameter (5T) that 
multiplies the calculated value of thrust using the above relations. 
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This parameter thereby acts as a throttle, increasing or decreasing the 
thrust as required for the ST bifurcation sweeps. 

Vehicle Characteristics. As with the majority of previous studies 
conducted on this subject, the vehicle geometry and aerodynamic 
characteristics used here will be basically the same as those used by 
Etkin (4:783-784). This allows nearly direct comparison of results which 
is helpful to check consistency of linearization at starting points and 
more importantly will be used to highlight the advantages and simplicity 
of using bifurcation analysis for even the simplest problems in 
atmospheric flight mechanics/dynamics. Etkin obtained his data from 
"simple Newtonian impact theory for a slender body (cone or wedge of 
about 3° semiangle) at moderate angle (4:784)." In this study a small 
change has been made to allow for lift at zero angle of attack, which is 
more representative of a hypersonic vehicle, however as seen in Figure 3 
Etkin's basic lift to drag ratio was followed fairly well. 

With these clarifications stated the geometry and aerodynamic 
characteristics are as follows: 


Geometry 

1 = 50 ft k = (I /m) 1/2 = 25 ft fca (I x -!,)/!,«-0.94 

W = 700,000 lb W7S = 30 psf (at sea level) 


Aerodynamics (dimensions are [rad'*] unless otherwise indicated) 



0.05 

0.50 

0.0133 

0.400 


'lO 

'»a 

;»q 

'i5bf 


0.00 

-0.0548 

-0.028 

-0.0822 [deg' 1 ] 
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C I / Cd 



Figure 3. Comparison of Cp/Cp with Etkin's work 
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II. Introduction to Bifurcation Analysis and Stability 


Equilibrium Solutions and Equilibrium Points 

At this point it is best to discuss some of the central points of 
bifurcation analysis to allow the uninitiated to understand what is being 
done, and to highlight what is being presented in the bifurcation diagrams 
and phase plane diagrams which follow. 

Continuation is the core on which the present application of 
bifurcation analysis is based. There are several types of continuation 
techniques, but discussing them is not appropriate here. What is required 
is a general description of what is obtained. The analysis here starts at 
an equilibrium point from which the continuation technique is begun. The 
equilibrium point is obtained by setting the time derivatives of the state 
variables to zero thus creating a set of homogeneous equations and solving 
simultaneously for the values of the state variables that satisfy the 
homogeneous equations. Once this starting point is obtained the 
continuation method can begin. The continuation method will search in the 
vicinity of this point until it finds another point which satisfies the 
set of homogeneous equations. This process continues over the given range 
of the specified parameter and within the bounds established for the 
states, until a complete parameter-dependent family of equilibrium 
solutions to the set of homogeneous equations has been found. An 
important feature of the bifurcation software AUTO is the ability to find 
the equilibrium solution path despite running into singular points (limit 
points and bifurcation points). AUTO uses a pseudo-arclength technique 
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which allows it to continue working despite encountering a singular 
Jacobian of the linearized set of equations. Other less robust techniques 
fail at these points where the slope of the solution with respect to the 
parameter is not unique, as in the case of a bifurcation, or undefined, as 
with a limit point. The reader is referred to the user's manual for AUTO 
and the text by Seydel for further information on this technique (3:12-16; 
12:116+). 

Simple Nonlinear Behavior: 

While knowing an equilibrium solution branch is interesting, the 
true value in this analysis for those interested in nonlinear effects is 
the accurate location of limit and bifurcation points. The reasons for 
interest in these singular points are many. In the case of a bifurcation 
this point represents the intersection of other solution branch(es). 
Depending on the type of bifurcating point there is the potential for 
complex motion arising from the nonlinear nature of the problem. In the 
analysis of longitudinal motion of a powered lifting hypersonic vehicle 
the two most prevalent singular points encountered are limit points, which 
may give rise to hysteresis type behavior or an exchange of stability, and 
Hopf bifurcations, which generally occur in this study when the phugoid 
mode eigenvalues cross the imaginary axis transversely and either lose or 
gain stability. Generally for the analysis of the longitudinal dynamics 
of a powered lifting hypersonic vehicle the Hopf point signals the loss of 
stability in the phugoid mode. A Hopf point is of special interest in the 
study of nonlinear dynamics as the behavior subsequent to a Hopf 
bifurcation is generally characterized by increasingly complex periodic 
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motion known as limit cycles. For systems with three or more degrees of 
freedom the Hopf bifurcation may also be the first step in the direction 
of chaos (12:25). 

While the mathematics of bifurcation analysis is not within the 
scope of this thesis an understanding of physical processes is. To 
further clarify some of the points made, and to provide a basis for 
understanding what a basic bifurcation study is all about, the following 
simple example is presented (10). 

Figure 4 shows a mass (m) connected by a mass-less rigid rod to a 
mass-less rigid sleeve which is around a spinning shaft. The friction 
coefficient between the spinning rod and the sleeve is constant therefore 
a constant torque is generated which is transmitted to the mass as a force 
(F) via the mass-less rigid rod. In-set in Figure 4 is the free-body 
diagram for the mass and the reference axes. Note the pendulum is held at 
a constant angular position by the torque applied to the sleeve. 

The scalar equation of motion for the mass m is: 

mib = F - mg sin(0) f 28 ) 


Where: F 

g 

m 

r 

T 


force generated by the constant torque (F=T/r) [lb] 

gravitational acceleration (ft/sec 2 ) 

mass (slugs) 

length of rod (ft) 

torque applied to sleeve (ft-lb) 


At a point of equilibrium the left hand side of eqn (28) is equal to zero, 
that is there is no change of the state variable when the forces acting on 
the mass are in equilibrium. 
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This brings out an important point to remember when looking at bifurcation 
diagrams: the stationary solution path is made up of equilibrium 
solutions (equilibrium points) and no change of the state variables is 
invol ved. 

For the equilibrium solutions eqn (28) is a homogeneous equation. 
In the case of the powered lifting hypersonic vehicle there is a set of 
homogeneous equations that are solved to locate the equilibrium solution. 
To actually conduct the bifurcation analysis a parameter must be 
established for which the continuation process finds a solution path. For 
the example problem, the parameter can be identified as A = F/(mg) and 
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the homogeneous equation can then be put in the form of eqn(29) below: 


0 e = sin' 1 (—) (29) 

* mg 

By varying A from -1 to 1, a plot of the equilibrium solutions for 6 e is 
easily obtained even with a hand calculator. The bifurcation analysis 
however, yields information on stability by linearizing the differential 
equation (or set of equations) and solving for the eigenvalues of the 
subsequent Jacobian thus providing the following bifurcation diagram. 

In looking at the bifurcation diagram it is seen there are at least 
two equilibrium solutions for each A in the open set (-1,1). The two 
points corresponding to A=-l or 1 are called limit points. It is clear to 
see that limit points have only one value of the equilibrium point for a 
given parameter and are points where the equilibrium solution path turns 
back with respect to the parameter; thus the slope is undefined. Note 
also that to one side of a limit point no equilibrium solutions exist yet 
on the other side two equilibrium solutions exist for each A. In the 
example problem the limit points also correspond to points where stability 
is either lost or gained, depending on the direction of the parameter A, 
however this is not always the case for limit points in general. Finally, 
note the convention used in bifurcation diagrams is to identify stable 
equilibrium branches with solid lines and unstable equilibrium branches 
with broken or dashed lines; other graphical conventions will be brought 
to the reader's attention as required. 
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Figure 5. Bifurcation Diagram for Example 
Problem 


The bifurcation diagram now shows part of its value in that the 
local, and at times global, behavior of trajectories can be predicted. 
For this example problem trajectories in the vicinity (domain of 
attraction) of a point on the stable equilibrium branches will oscillate 
about the equilibrium point since no damping is included. If damping were 
added the trajectories would be attracted to the equilibrium point. This 
can be visualized in thinking of the mass at an angle 9 g that corresponds 
with the stable branch for a given value of X. Given a small perturbation 
in any direction the mass will return to the equilibrium point, that is it 
will be attracted to the stable branch. Trajectories in the vicinity 
(domain of repulsion) of a point on one of the unstable solution branches 
will be repelled. Again think of the pendulum mass at an angle © e that 
corresponds to a point on one of the unstable branches. Given a small 
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perturbation the mass will not return to its original position or 
trajectory but will move away and in the example problem will swing around 
the shaft. Note finally, that the domain of attraction can be limited. 
In the example problem, even if the mass is in a stable position in the 
lower part of the pendulums arc, if given sufficient perturbation in the 
direction of the force F the mass can be made to swing around the shaft 
and in the absence of any damping action will not settle back to a stable 
position; i.e. the domain of attraction is limited. Finally, note for 
values of X > 1 and A < -1 the pendulum will swing continuously about the 
shaft - - no equilibrium exists. 

The foregoing example has provided a good picture of the very basics 
of bifurcations analysis and some of the phenomena that may occur when 
analyzing nonlinear problems. In addition it has provided some physical 
significance to limit points. It was however limited in its ability to 
fully demonstrate possible behavior as it has only one degree of freedom. 
Problems with three or more degrees of freedom can develop many other 
phenomena and a variety of ways to exchange stability. 

Limit Cycles (Orbits) 

Hopf bifurcations, which are a central feature in the study at hand, 
give rise to periodic solution branches which on a bifurcation diagram are 
shown by plotting the maximum and/or minimum values of limit cycles 
(Figure 6). These surfaces projected in the phase plane, or phase space 
for systems with degrees of freedom greater than two, are called limit 
cycles or orbits and represent a surface of periodic solutions surrounding 
a equilibrium solution branch (Figure 7). 
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Figure 6. Limit of Periodic Motion about a 
Solution Branch 
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Just as with equilibrium solution branches, periodic solution 
branches can be stable or unstable and will either attract or repel 
trajectories within their domain of influence (12:12-25). Local stability 
of the periodic solutions is based on Floquet theory involving the 
monodromy matrix (12:240-248). For the purpose of this study one need 
only know the relationship of the eigenvalues of the monodromy matrix with 
the criteria for local stability of the periodic solution branch. The 
monodromy matrix always has one eigenvalue (Floquet multiplier) in the 
complex plane at z=l. This is subsequently used as a test for accuracy in 
calculating the other multipliers in AUTO. Through establishing a 
relationship between the remaining Floquet multipliers and a Poincard map 
or return map it is determined that if the modulus of the remaining 
Floquet multipliers are each less than unity then the periodic solution is 
stable (i.e. attracting). AUTO computes the Floquet multipliers and uses 
this to determine stability. As mentioned above the accuracy of this 
calculation must be checked based on the one Floquet multiplier at 2=1. 
Since five ordinary differential equations make up the full set of 
equations needed to model the longitudinal dynamics for this study then 
for the periodic branches to be stable four Floquet multipliers must each 
have a modulus less than unity and one must have a modulus equal to unity. 
Further discussion of nonlinear behavior is beyond that required to begin 
the analysis, these concepts will be expanded on as the results of the 
study are examined. 
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Two Parameter Continuation: 


Upon locating a Hopf bifurcation point or limit point AUTO can be 
used to perform a two parameter continuation which will show the evolution 
of bifurcation points as the second parameter is varied; ordinary 
equilibria is ignored. For the study here the two parameter continuation 
will be used to determine the value of a gain in a feedback loop to 
attempt to stabilize a system. Specifically by plotting one parameter 
against another a curve of limit points or Hopf points is generated. If 
for instance, the first parameter is the body flap deflection or a pitch 
command and the second parameter is the pitch rate feedback gain in a 
pitch attitude feedback loop, the plot of the two parameters would show 
how the Hopf point (the point of stability exchange) moves with the pitch 
rate gain. A typical plot might look something like Figure 8. 

In this figure one can see that for values of pitch rate gain the Hopf 
point will move along the curve. If one wishes to delay the occurrence of 
the Hopf point relative to the magnitude of the body flap deflection 
(thereby delaying the exchange of stability) one would raise or lower the 
gain as appropriate until the Hopf point is moved far enough to meet the 
system requirements based on the body flap deflection. Some simple 
feedback techniques are examined further in Section IV. 
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Figure 8 


Curve of Hopf Points for Two Parameter 
Continuation with Pitch Rate Gain and 
Body Flap Deflection 
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III. Bifurcation Analysis of Longitudinal Dynamics 


To begin the actual analysis the set of ordinary differential equations 
given by eqns(7,9,10,13 and 16) must be solved for an equilibrium point. 
As stated before this equilibrium starting solution is required for the 
continuation method to begin. Taking the equations of motion and 
substituting the force and moment coefficients yields: 


T , v PV 2 SC D . , . 

V = — cos(a) - —- - gsin(y) 

m 2m 


(30) 


4 • 1 - -~sin(a) - ^ + SM-cosiy) 

mV 2m i 


(31) 


5 = g + —cos (y) 
r 


(32) 


= P^ . c . - |J[[i£li£]sin(20) 

2 ly 2 *{ J y ) 


(33) 


f = Vsin(y) 


(7) 


where the force and moment coefficients are as defined before in eqns(20, 
21 and 22) and are again repeated below for convenience. 


+ C x j, a 2 


( 20 ) 


~ CL0 + Q* 0 


( 21 ) 
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CjB = ^mo + ® + ^XXJ + 


(22) 


The equilibrium values are obtained by first setting the left hand side of 
equations (7, 30, 31, 32 and 33) to zero and solving for the states (V, a, 
q, 6, and r). To make things easy y is set to zero (y = 0), which means 
a = 9 = 0 . From this point the following equations are obtained that 
describe the requirements for equilibrium with r = r c , a = 9 = 0: 


Yo = 0-0 


(34) 




[ 9 (r 0 )r 0 

1+ 

fp(r 0 )ar 0 C^ 

l 2 m j 


So = - 


Yi 
10 


(35) 


(36) 


T(r 0 ) = -| p {r 0 )SC D vt 


(37) 


C m = 0.0 

in 


(38) 


This set of equations is programmed into the user provided subroutine FUNC 
which AUTO requires to find the flow of equilibrium solutions as a 
specified parameter is varied. Note the body flap parameter (5bf) shows 
up in the last term of eqn(22). Also of interest is the value of the 
thrust at the starting equilibrium point - - it equals the drag as was 
discussed earlier for the thrust laws. The AUTO software package was used 
to develop the equilibrium solution branches and identify bifurcation 
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points for the system described by eqns(7, 30, 31, 32 and 33). The 
parameters used to conduct the analysis were a body flap parameter (6bf) 
and a parameter which scaled the thrust (6T); simulating a throttle. This 
study was conducted for several starting equilibrium altitudes, which 
amounts to selecting a starting altitude (r = rj), and setting the angle 
of attack and pitch angle to zero and letting AUTO solve for the velocity, 
pitch rate and thrust. The starting altitudes ranged from 50,000 ft to 
700,000 ft. Again note that for each starting equilibrium altitude a 
different value of the initial equilibrium thrust is obtained since, as 
discussed before, thrust must equal drag at the starting equilibrium 
point. Note that the thrust decreases with increasing altitude. The 
resulting bifurcation diagrams are identified by their starting altitude 
as the thrust level, which is fixed by the drag at the starting altitude, 
makes a difference in the behavior of the system. 

One can see the equations of motion are clearly nonlinear with the 
states all interrelated, However two states exert the primary influence 
by virtue of the type of problem being analyzed, these are the velocity 
and the radius. It is the way these two states vary to achieve 
equilibrium and the fact that the behavior being observed is the behavior 
of the equilibrium solution path that makes for results that are not 
intuitively obvious relative to the effects of changing the body flap or 
throttle (or more precisely thrust variations). This point must be 
emphasized since most traditional dynamists are used to dealing with the 
time history of trajectories or frequency response given some control 
input. As discussed in section II the bifurcation diagrams which are used 
here are not time histories but a collection of equilibrium points that 
provide the value of the states relative to a parameter. 
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The fact that the solution branch is made up of equilibrium points implies 
that no change of the state variables occurs at the individual points 
along the curve; this can at times cause confusion so beware. 

Three cases were initially considered based on the three separate 
thrust laws discussed in previous sections. The results of the rocket 
whose thrust varied with altitude differed little from the more standard 
constant thrust rocket used by most previous authors working on this topic 
so the case was dropped. The remaining two cases were studied to 
investigate the behavior the rather simple nonlinear model would generate. 
What follows is a discussion of the results. Since the effect of the body 
flap and throttle are significantly different it is best to discuss them 
separately. 

Body Flap Parameter (5bf) Variation 

The body flap parameter mathematically represents a deflection of 
the body flap in degrees. This value is changed to radians and affects 
the value of C N via eqn{22). A change in the value of C N generally causes 
a change in the vehicles angle of attack and pitch attitude. It is 
interesting to note that the primary influence of the body flap on the 
behavior of the equilibrium solution path is, that in changing the angle 
of attack the value of the lift coefficient is changed which is a key 
parameter in establishing the value of the velocity for a equilibrium 
orbit (16:321-344). Therefore it is best to think of the body flap as a 
control by which lift is modulated. 

To understand the equilibrium solutions obtained relative to the 
body flap parameter one must examine the relationship of the velocity and 
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radius with C- b . The following equation relates the velocity to and 
radius and is central in understanding what is occurring when looking at 
a equilibrium solution path displayed in a bifurcation diagram. 




N 


rig - 
1 + 


rsin(oc) 

m 

p5Cj,r 

2m 


) 


(39) 


This equation results from the homogenous set of equations for equilibrium 
and represents the velocity - altitude relationship for everything from an 
unpowered satellite or lifting vehicle in equilibrium orbit to the case 
here of a powered lifting vehicle in equilibrium orbit. The velocity 
plays a key role in providing forces sufficient to balance the weight of 
the vehicle. At high altitudes the centrifugal force is the primary 
means by which the vehicle balances the weight, where at lower altitudes 
the lift, which is a function of velocity is the primary balance to the 
weight. The way m which the equilibrium solution path moves in order to 
balance all forces is quite interesting and as stated before not always 
obvious. Those interested in knowing more about how velocity, altitude 
and the lift coefficient are related in orbital flight are refered to 
reference ( 16 ). 

The investigation was conducted by performing a continuation from 
the equilibrium starting point (i.e. r=r 5 ) using the body flap parameter 
to control the initial direction of the continuation. A body flap sweep 
is defined as the summation of the equilibrium branches obtained from 
performing the continuation in the directions associated with a positive 
flap deflection and a negative flap deflection (control surface movement 
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downward being positive). The flap was constrained to +/- 90 degrees by 
fixing the limiting values for the parameter in AUTO. Starting altitudes 
for the constant thrust rocket case were varied between 100,000 ft and 
700,000 ft and for the air-breathing engine case ranged from 50,000 ft to 
400,000 ft. Note in all cases as the starting altitude is increased the 
starting equilibrium value of the thrust is decreased. The equilibrium 
solution paths were tabulated, stability determined and all simple 
bifurcations, limit points and Hopf bifurcation points were located. 
Having mapped out the equilibrium solution and identified the various 
singular or bifurcation points any Hopf bifurcations were continued to 
obtain the limit cycles. This process is accomplished by taking the 
equilibrium conditions at the point identified as a Hopf bifurcation as a 
starting point for AUTO's continuation method. Specific software routines 

t 

in AUTO are used to perform the required functions to obtain the periodic 
solution branch. These data are generally interpreted graphically to 
obtain a general feel for the local behavior of trajectories in the 
vicinity of the solution branches; these graphs are known as bifurcation 
diagrams. Since bifurcation diagrams are meant to convey information 
about the behavior of the system, it seems only natural that a method or 
convention be established for presenting data on these. The reader is 
encouraged to take note of the following rules for conveying information 
about the types of solution branches and their local stability 
characteristics. 

Equilibrium solution branches are presented as lines. Solid lines 
indicate stable solution branches and any type of broken or dashed 
lines indicate an unstable equilibrium solution branch. 

Periodic solution branches are shown as circles or dots which 
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generally indicate the maximum amplitude of the periodic motion. 
Stable periodic solution branches are shown as solid dots or filled 
circles. Unstable periodic solution branches are indicated by open 
circles. 

These conventions are adhered to throughout this paper. 

Constant Thrust Rocket Case. The following bifurcation diagrams 

were generated using AUTO as described previously, however further 

explanation of the way the constant thrust value is determined. As stated 

several times before, and shown explicitly in eqn(37), the thrust equals 

the drag at the starting equilibrium solution with y = 0 radians. For 

the constant thrust rocket case this value obviously is fixed over the 

entire body flap sweep. This means that two equilibrium solution points 

with the same altitude but obtained from body flap sweeps starting at two 

different altitudes will not have the same value of thrust and in general 

will have different values for the other states as well. Note finally, as 

starting altitudes are increased the value of the thrust decreases. 

Figures 9 and 10 show a collection of bifurcation diagrams for each 

state (note a = 0 ) for the body flap (5„<) sweeps from 100,000 and 300,000 

•» * 

ft respectively. While the behavior is nonlinear there is not a great 
deal of interest occurring over most of the equilibrium branches. 

Figure 9 is characteristic of the constant thrust rocket case with 5^ as 
the parameter for starting altitudes less than 150,000 ft and Figure 10 is 
characteristic of the constant thrust rocket case with 5^ as the parameter 
for starting altitudes between 150,000 ft and about 360,000 ft. One of 
the first things to note in each figure is the system is unstable 
(indicated by the dashed line). 
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Figure 9. Bifurcation Diagrams for Body Flap Sweep from 100 kft 
Constant Thrust Rocket 
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Figure 10. Bifurcation Diagrams for Body Flap Sweep from 300 kft 
Constant Thrust Rocket 










The instability is caused by the nonoscillatory height mode and was the 
case for both thrust laws representing a rocket. This behavior is 
consistent with the previous studies and follows as a result of the way 
the sum of the longitudinal forces change with respect to altitude 
(13:472; 7:283). Interestingly enough the height mode would generally 
stabilize at some point for starting altitudes between about 360,000 ft 
and 530,000 ft and would occur associated with a limit point. 
Frustratingly this is about the point where the mode normally thought of 
as the phugoid mode (based on the longer period) would go unstable. 

From the stand point of nonlinear analysis not much of interest is 
occurring for the body flap sweeps with starting altitudes less than about 
360,000 ft. Specifically, the system is unstable with generally no limit 
points from which jump phenomena may occur and no simple bifurcations. 
Only if the continuation of a equilibrium branch associated with 
increasing altitude is allowed to go long enough, to where the aerodynamic 
pitch damping is lost due to the very low density at very high altitudes, 
will a Hopf point be found. For the body flap sweep from 300,000 ft 
(Figure 10) the Hopf point occurs at 5^ = +/- 57.9° , and an altitude of 
615,120 ft for the branch associated with a negative body flap deflection 
and 614,970 ft for the other branch; this symmetric behavior is not as 
closely followed at lower starting altitudes where the thrust levels and 
densities are greater. 

Concentrating now on the unstable periodic branches found from the 
body flap sweep from 300,000 ft one can see from the bifurcation diagrams 
in Figure 10 that as the unstable periodic branch progresses the amplitude 
of the limit cycles become fairly substantial, however the period is on 
the order of 5000 seconds and since the periodic branch is unstable 
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trajectories would not approach it. Recall that the limit cycles of a 
periodic branch, indicated by circles, show the maximum and/or minimum 
amplitude of the ■; * '.die motion. The stability of the periodic branches 
is interesting in that it is just barely unstable with a conjugate pair of 
Floquet multipliers just outside the unit circle. The height mode at the 
point where the phugoid mode goes unstable is very near the imaginary axis 
with a time constant on the order of 10^ seconds; this is consistent with 
previous studies (4:786). 

Figure 11 shows the time history over one period for the right 
periodic branch with 5^- = 56.73° and 5^= 42.25°. Figure 12 shows the 
time history over one period for the left periodic branch with 5 hf = -56.6° 
and 6„ ; = -41.8°. Note that near the Hopf point on either branch the 
motion is slight (ie just leaving equilibrium) but as the parameter is 
changed to move along the periodic branch the motion increases. Once 
again the behavior of the periodic branches shown in this body flap sweep 
(from 300,000 ft) is characteristic of the behavior of the periodic 
branches occurring from body flap sweeps starting at "lower" altitudes 
that subsequently extend to altitudes above 500,000 ft where pitch 
stability is lost. 

In order to provide a complete look at the periodic behavior of the 
limit cycles, as well as provide a connection with more classic 
longitudinal analysis, the limit cycles are projected into the phase plane 
with the flight path angle. Figure 13 shows two limit cycles from the 
right periodic branch; one for 6^ = 56.7° and one for 5^ = 42.25°. 

Figure 14 shows two limit cycles from the left periodic branch; one at 
6.,- = -56.64° and one for 6^ = -41.8°. In looking at these figures one 
sees clearly that motion is associated with the periodic branches. 
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Figure 11. Left Branch Limit Cycles for 6^ = -56.6°and -41.8° 
Body Flap Sweep from 300 kft : Const Thrust Rocket 





Figure 12. Right Branch Limit Cycles for 6^ = 56.73° and 42.25° 
Body Flap Sweep for 300 Kft : Const Thrust Rocket 
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Figure 13. 


Right Branch Limit Cycles in Phase Plane for = 56.73° 
and 42.25°. Sweep from 300 kft : Const Thrust*Rocket 
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Figure 14. Left Branch Limit Cycles in Phase Plane for 5^ = -56.64° 
and -41.8°. Sweep from 300 kft : Const Thrust Rocket 
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Also, the behavior of the translational states is very much like the 
classic phugoid mode, i.e. showing a steady exchange of potential and 
kinetic energy. While the periodic branches are unstable, the growth of 
the nonlinear behavior in velocity and altitude remains sinusoidal.In 
terms of the rotational states, one sees somewhat more complex behavior 
which should be expected since the periodic branches arose from the loss 
of pitch stability at the very high altitude. 

The most interesting body flap sweep for the constant thrust rocket 
case resulted from the starting altitude of 400,000 ft. The body flap 
sweep bifurcation diagrams and expanded views for a and altitude are shown 
in Figures 15, 16 and 17. The behavior of the limit cycles of the 

periodic solution branches are certainly visually interesting. Note the 
periodic branches contain several limit points which explains their 
complex twisting about. 

The Hopf bifurcation occurred at 5 b£ " +/- 8.93° and generally 

speaking is not of great significance since the equilibrium branch was 
unstable to begin with and the periodic branch starts out unstable and 
encircles an unstable equilibrium branch. On closer inspection of the 
periodic branches one will see (Figures 16 and 17) that there is a portion 
of both periodic branches, starting at = +/- 7.82° and continuing to 

5 bf = +/, “ 6 - 05 °' that gains stability by the crossing of a conjugate pair 

of Floquet multipliers. This type of stability exchange is associated 
with bifurcation to a torus. What this implies is that trajectories 
within the domain of attraction of the periodic branch will be drawn into 
periodic motion with two frequencies; one describing the component of the 
motion in the circumference direction of the torus and one around the 
cross section of the torus. The path of the trajectory can be visualized 
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Figure 15. Bifurcation Diagrams for Body Flap Sweep from 400 kft 
Constant Thrust Rocket 


as spiralling around the inside of an inner tube (i.e. torus) in 
hyperspace (12:263,264). 

The accuracy of the solution is very good with the Floquet 
multiplier that is supposed to be equal to one (z=l), precisely equal to 
one. The fact that stability is gained then lost so quickly indicates 
that one or more Floquet multiplier(s) is(are) very near the edge of the 
unit circle and upon inspection of the output from AUTO one finds this to 
be the case with one pair inside the unit with modulus = 0.99884 and a 
second pair with modulus = 0.98. This would indicate a weakly attracting 
limit cycle for the range of 5^ where the periodic branches are stable. 
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6bf Sweep from 400 Kft : Const Thrust Rocket 













5M'--Sweep from 400 Kft : Const Thrust Rocket 



Figure 17. h Bifurcation Diagram for Body Flap Sweep from 400 kft 
Constant Thrust Rocket 











Figure 18 shows the motion of the limit cycle over one period for 
periodic branch 1 (left branch) for a point in the stable region at 
Sjjj = -6.83° and a point near the end of the calculated portion of this 
periodic branch at 5^ = -2.89°. Figure 19 shows the motion of the same 
two limit cycles (5^ = -6.83° and 5 bf - -2 .89°) over one period in the 
phase plane with the flight path angle. Figure 20 shows the motion of the 
limit cycle over one period for periodic branch 2 (right branch) with 
5 b £ = 8.73° in the stable region and for 5 bb = 3.57° near the end of the 
calculated branch. Figure 21 shows the limit cycles at 6 bb = 8.73° and 
5^ = 3.57° in the phase plane with flight path angle. Observing the 
behavior of the limit cycle over one period at several points like this 
shows why the nonlinear behavior associated with period solutions are so 
interesting. Looking at Figures 18 and 20 one can see a kind of wave 
changing in amplitude as the parameter is varied. Figure 22 shows 
qualitatively the growth of the nonlinear behavior in a as periodic 
branch 2 grows. It is this type of behavior, for systems with three or 
more degrees of freedom, that leads to more fascinating subjects like 
chaos and the Hopf point is as Seydel puts it, "the door which opens from 
the small room of equilibria to the large hall of periodic solutions 
(12:61)." 

On a somewhat different note, the altitude where the Hopf point 
occurs or. both branches is just about 450,000 ft. This is in the altitude 
range where Etkin determined, and later others modified, that the period 
of phugoid and pitching modes came very close to each other (4:787-788; 
15:17-20). The general conclusion of these earlier works is that there 
would be significant coupling of the two modes at this so called 
"resonance altitude" (15;7). 
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18. Limit Cycles for Periodic Branch 1 (5^ =-6.826° & -2.89°) 
Body Flap Sweep from 400 kft : Const Thrust Rocket 
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Figure 19. Limit Cycles in Phase Plane, Periodic Branch 1 

<6l£ = -6.826° & -2.89°) : Constant Thrust Rocket 
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Figure 20. Limit Cycles, Periodic Branch 2,(5^ = 8.738° & 3.57°) 
Body Flap Sweep from 400 kft : Const Thrust Rocket 






Figure 21. Limit Cycle in Phase Plane, Periodic Branch 2 (5^ = 8.738° 
& 3.57°), Sweep from 400 kft : Const Thrust Rocxet 
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Figure 22. Growth of Nonlinear Waveform, Periodic Branch 2, 

a Bifurcation Diagram, Body Flap Sweep from 400 kft 


In looking at the limit cycles in the preceding figures there does 
not seem to be the strong coupling predicted. Some coupled motion is 
evident in that the rotational states (a, 0, and q) go through sub¬ 
oscillations in each overall period while the translational states 
maintain a sinusoidal motion with a very regular period. Since the period 
is nearly that of the circular orbit for the same geocentric radius, it 
would seem that what is observed is a barely unstable elliptical orbit 
with the vehicle pitching about its y-axis at some rub-frequency greater 
than the frequency of the orbit (i.e. overall frequency of the limit cycle 
for the given parameter). 

As a final note, notice the limit points on the equilibrium branch 
where the equilibrium solution path changes direction. 




It is interesting in that it exists and is associated with the height mode 
stabilizing. Notice how the direction of the body flap deflection changes 
in order to maintain the original direction of the equilibrium solution 
path. 


Air-breathing Engine Case. The procedure for continuation and 
subsequent analysis for the air-breathing case with the body flap as the 
parameter, was the same as that for the constant thrust rocket analysis. 
The results are somewhat more interesting in that the height mode is 
generally stable for the air-breathing case below approximately 
380,000 ft thus the entire system is stable at these "lower" altitudes. 

A phenomenon of little physical significance, but interesting 
nonetheless for the air-breathing case with body flap sweeps starting from 
equilibrium points between approximately 100,000 ft up to approximately 
360,000 ft is that the velocity goes very nearly to zero for the portion 
of the body flap sweep that has a negative body flap deflection. Before 
the velocity actually gets to zero, a Hopf bifurcation occurs, then a 
simple bifurcation which has two solution branches. Of the two branches 
one stable and back-track the original equilibrium solution branch for all 
the states, and the other is unstable. The unstable branching solution 
back-tracks a, and altitude, but takes the negative of it's original value 
for velocity and pitch rate. Figure 23 shows this bifurcation diagram for 
the equilibrium branches only. Further work is needed to finish exploring 
this behavior. Notice that in Figure 23 the + - symbol indicates the 
stable branching solution and the x - symbol indicates the unstable 
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branching solution. Finally the original branch also turns back on itself 
after going to nearly zero subsequent to the Hopf point; this back¬ 
tracking original branch turns back as an unstable branch but regains 
stability as it repasses the Hopf point (as expected). 






Figure 23. Bifurcation Diagram for Body Flap Sweep from 100 kft 
Air-Breathing Engine 


For ease of discussion, Figure 24 shows all of the positive body flap 
sweep from 100,000 ft, but only the stable portion of the negative body 
flap sweep from 100,000 ft for the air-breathing engine case. This 
diagram contains the Hopf bifurcation but none of the "back-tracking" 
solutions. A big difference from the previous case is readily apparent, 
the system is stable up to the Hopf point at which time the phugoid mode 
loses stability and note how low the altitude is ( approximately 73,000 
ft). Notice this is a subcritical bifurcation (12:72); that is an 
unstable periodic branch encircles a stable equilibrium branch. What this 
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implies is trajectories near the periodic solution branch but within the 
domain of attraction of the equilibrium branch will be drawn to the 
stable equilibrium branch and generally no further changes of the states 
will occur. 




Figure 24. Partial Body Flap Sweep from 100 kft for the Air-Breathing 
engine case 


Figure 25 shows the limit cycles over one period for a point near 
the Hopf bifurcation {5^ = -52.77° ; T = 77 sec) and the point on the end 
of the calculated periodic branch (5 b{ = -23.54° ; T = 140 sec). Clearly 
visible is the increase in nonlinear behavior as one moves along the 
periodic branch. Near the Hopf point one can see the motion is nearly 
constant and what little variation is there is sinusoidal. For points 
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Figure 25. Limit Cycles for = -52.77° & 6^ = -20.54°, Body Flap 
Sweep from 100 kft : Air-Breathing Engine 


farther along the periodic branch the nonlinear behavior, due to the 
underlying nonlinear equations, truly begins to blossom. 

Figure 26 shows the limit cycles for 5^ = -52.77° and -23.54° in the 
phase plane. Once again the classic phugoid-like behavior of the 
translational states is seen. Note as before in the constant thrust case 
the rotational states have this sub-oscillation. In contrast to what was 
seen in the constant thrust rocket case, here the phugoid mode has gone 
unstable at relatively low altitude (approximately 73,000 ft). 

A body flap sweep from 400,000 ft for the air-breathing engine case 
is shown in Figure 27. The general characteristics for this sweep are the 
same as for the constant thrust rocket. However here the periodic branch 
remains unstable. As with the constant thrust rocket case, each 
equilibrium branch has a limit point where the height mode stabilizes. 
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night Poth Anglo [dag] night Path Angla [dag] 

Figure 26. Limit Cycles in Phase Plane (5^ = -52.77° & -23.54°) 
Body Flap Sweep from 100 kft : Air-Breathing Engine 


There are also several limit points found along each periodic branch; just 
as with the constant thrust rocket body flap sweep from 400,000 ft. 

Figures 28 and 29 show the limit cycles over one period and the 
limit cycles in the phase plane for the right periodic branch with 
5 fc£ = -0.00078° and 0.0118°. Figures 30 and 31 show the limit cycles over 
one period and the limit cycles in the phase plane for the left periodic 
branch with = -0.0046° and -0.011°. These points are as before, used 
to show the behavior of the limit cycle near the Hopf point versus the 
behavior farther down the periodic branch. The same basic behavior is 
seen for the translational states as discussed for previous cases. 
However of interest is the phase shift occurring along the left periodic 
branch. It is interesting to see how the 3°ft branch differs markedly 
from the right branch even though the two look relatively symmetric. 


57 












Body Flap Deflection [deg] Body Flap Deflection [deg] 


Figure 27. Collection of Bifurcation Diagram for Body Flap Sweep from 
400 kft : Air-Breathing Engine Case 
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Figure 28. Limit Cycles over one Period 5bf= -0.0007° and 0.0118° 

Body Flap sweep from 400 Xft : Air-Breathing Engine Case 





Figure 29. Limit Cycles in Phase Plane for 6bf= -0.0007°and 0.0118° 

Body Flap Sweep from 400 kft : Air-Breathing Engine 
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Figure 30. Limit Cycles, Left Branch,( 6^ = -0.0046°& -O.Olf) 
Body Flap Sweep from 400 kft : Air-Breathing Case 





Figure 31. Limit Cycles in Phase Plane, Left Branch (5^ = -0.0046° 
& -0.011°), Sweep from 400 kft : Air-Breathing Engine 
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Throttle Parameter (5T) Variation 


The throttle parameter acts as a multiplier to scale the equilibrium 
value of thrust. In the case of the constant thrust rocket it scales the 
value of the thrust at the starting altitude which was set equal to the 
drag. For the air-breathing engine case the throttle parameter scales the 
thrust as with the rocket case (i.e. which was set equal to drag at the 
starting altitude), however for the air-breathing case the thrust varies 
as a function of altitude via T = T 0 [p(r)/pjj]. A change in thrust causes 
a corresponding change in drag to maintain equilibrium so once again the 
velocity and radius begin to play an important role. However an 
interesting characteristic of the equilibrium solution paths where the 
throttle is the parameter is the lack of any modulation in the coefficient 
of lift or drag; which remain effectively constant at all throttle 
settings and radii. 

The analysis using the throttle as a parameter was done in the very 
same manner as the body flap parameter. As before a starting altitude was 
selected (with a = 9 = 0 radians) from which the equilibrium values for 
the remaining states were calculated. From this point the throttle 
parameter was first increased from 6T=1.0 then decreased. This yielded 
a complete throttle sweep. 

Constant Thrust Rocket Case. Not much happened with this case. From 
Figures 32 and 33 one can see the system is unstable, due to the height 
mode, and is nonlinear but no bifurcations were detected. It appears 
without the pitching associated with the body flap there is little to 
drive the nonlinear nature of the problem. In contrast to the cases where 
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the body flap was used as the parameter the angle of attack in 
cases is basically zero. 
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Figure 32. Collection of Bifurcation Diagrams for Throttle Sweep from 
100 kft : Constant Thrust Rocket 
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Figure 33. Collection of Bifurcation Diagrams for Throttle Sweep from 
400 kft : Constant Thrust Rocket 
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Air-Breathing Engine Case . For the air-breathing case the most 
interesting behavior was at altitudes below 400,000 ft. The following 
discussion characterizes the general nature of what was found at these 
"lower" altitudes. Figure 34 shows the 5T sweep for the air-breathing 
engine from 100,000 ft. For the portion of the sweep where 5T < 1.0 the 
system remains stable. For the portion of the sweep where 5T > 1.0 the 
system loses stability subcritically at a limit point (5T = 18.99 ) where 
the height mode crosses the imaginary axis. Just after this occurs (5T = 
18.93) a Hopf Bifurcation point is found. Since the limit point preceded 
the Hopf point this bifurcation is not classified either supercritical or 
subcritical and as before is really of little physical value other than to 
perhaps give an idea of the bound on the allowable perturbation to remain 
in the vicinity of the equilibrium solution branch. 

Figure 35 shows an expanded view of the area around the Hopf point 
and Figure 36 shows just the maximum limit cycles from the bifurcation 
diagram. These limit cycles show quite a variation of amplitude for the 
rotational states along the periodic branch as 6T is increased. Looking 
at Figure 37 one sees smooth sinusoidal behavior in the translational 
states even though the limit cycles are for points well toward the end of 
the calculated portion of the periodic branch (5T= 19.063 and 19.078). 
The rotational states however, display some relatively high frequency sub¬ 
oscillations. Although the amplitudes of these sub-oscillations are quite 
small and the period is on the order of 150 seconds (see Figure 38). In 
examining the phase plane representations of the limit cycles versus the 
flight path angle (Figures 39 and 40) one sees the translational states 
clearly displaying the motion that can be associated with an elliptical 
orbit. Further support of this view is given from Figure 41, where the 
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variation in the overall parameter dependent period of the periodic branch 
is shown relative to the circular orbital period for the values of the 
states at the given 5T. 
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Figure 34. Collection of Bifurcation Diagrams for Throttle Sweep from 
100 kft for Air-Breathing Engine Case 


64 






I’lli li Hull- | Dvk/SitJ Velocity (Ft/Sec] 



Thrust Scaling Factor Thrust Scaling Factor 


Figure 35. Expanded View of Region near Hopf Point for 5T sweep 
from 100 kft for Air-Breathing Engine Case 
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Figure 36. Maximum Limit Cycles for the 5T sweep from 100 kft for 
Air-Breathing Engine Case 
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38. Expanded View of a Limit Cycle, 5T = 19.063 Throttle Sweep 
from 100 kft : Air-Breathing Engine 
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Figure 39. Limit Cycles in Phase Plane for 5T = 19.063 
Air-Breathing Engine Case from 100 Kft 
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Figure 40. Limit Cycles in Phase Plane for 5T = 19.078 
Air-Breathing Engine Case from 100 kft 



Figure 41. Comparison of Circular Orbital Period and Period of Limit 
Cycles for Throttle Sweep from 100 kft : Air-Breathing 
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IV. Model Stabilization 


Simple Feedback Options 

In Stengel's work (13) he presented a summary of several possible 
feedback schemes using several controls. Table 1 in section I is a 
reproduction of the table found in Stengel's paper. In looking at these 
possible schemes one becomes aware of the dichotomy regarding the height 
mode and the phugoid mode. Most stabilizing feedback for the height mode 
destabilizes the phugoid and vice versa. In Berry’s work (2) similar 
feedback options were tried for the linear approximation technique and 
found to display the same behavior as found by Stengel. This "inverse" 
relationship between the height mode and the phugoid mode plus the 
restrictions imposed by the simplicity of the vehicle model made it beyond 
the scope of this thesis to actually stabilize the height mode. The 
success in stabilizing the height mode and not destabilizing the phugoid 
lay in developing a control law/technique to properly modify the way the 
longitudinal forces vary with height and velocity (as discussed in 
sec.ion I). Minor success was experienced m dealing with the phugoid by 
using pitch rate feedback to the body flap. It should be noted that pitch 
rate feedback in general has little afuc- on the phugoid roots however it 
was a technique that could be easily managed and did demonstrate the 
concept. For a given change in the pitch rate feedback the phugoid mode 
could be improved, but so slightly that in a practical sense it was 
worthless. Figure 42 shows the basic feedback loop with the pitch rate 
relative to the earth fed back in a negative feedback loop to the body 
flap (5 b{ ). 
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Figure 42. Pitch Rate Feedback Loop for Model Stabilization 


Note the parameter now becomes the value input as the command (5^). 

The new value of 5^ for use in eqn(22), the moment coefficient equation 
is given by the following: 

(«- d) ^ <»> 

where: 5^ = new parameter for body flap sweep control (deg) 

Kq * pitch rate feedback gain 

Figure 43 is sufficiently representative of all the cases where pitch rate 
feedback was used. It must be emphasized that the curve in Figure 41 is 
a curve of Hopf bifurcation points. What can be seen is that as the value 
of Kq is changed the location of the Hopf point relative to the body flap 
deflection is changed. In this case a gain of Kq=20 pushes the Hopf point 
down the curve the farthest. However, as one can clearly see the 
improvement is extremely small; to the point of being basically 
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constant thrust rocket system and attempting to use alternative feedback 
of other states, it seems that without recasting the mathematical model to 
allow for more reasonable feedback control, say with attitude command 
inputs, significant stabilizing routines cannot be obtained. 



Figure 43. Movement of Hopf Bifurcation given Pitch Rate Feedback to 
the body flap. Air-Breathing Engine from 100 kft 
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V. Conclusions and Recommendations 


Conclusions 


Overall this study has not produced results that would be considered 
"Earth shaking" or significantly different from previous work. It has 
however resulted in some worthwhile accomplishments, not the least of 
which is demonstrating how easy it is to use bifurcation analysis on 
problems related to hypersonic flight. This method provided much the same 
information as obtained by other authors using perturbation techniques, 
yet gave a much greater view of the actual effects of the nonlinear nature 
of the problem. 

In terms of comparison to previous work, it was found that the 
period and damping of the phugoid and pitching modes was similar to the 
behavior discovered by Vihn and Dobrzelecki and verified by Markopoulos, 
et al. (15:16-18; 7:286,287). Their study showed that the two modes do 
not cross for the linearized model as Etkin had concluded, but instead 
come very close together then diverge (4:785,786); this was the case here 
as well. From past work the behavior associated with the two modes in the 
vicinity of the "resonance altitude" (15;16) are expected to be coupled 
and behave nothing like that expected cf the classic pitching and phugoid 
modes. It was shown here that there is significant departure from the 
classic behavior of the rotational states, in that they show significant 
motion when the phugoid mode goes unstable. However the translational 
states act basically as expected. Looking at the limit cycles of the 
periodic branch associated with the "resonance altitude" one sees what 
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could be described as a one way coupling from the translational states to 
the rotational states, with the rotational states experiencing 
sub-oscillating for each period of the translational states. In addition 
the loss of stability at high altitude which Etkin concluded would occur 
due to the loss of pitch damping and a destabilizing moment due to the 
effect of the gravity gradient was seen (4:785,787). 

For the body flap bifurcation analysis the most interesting findings 
are the results associated with the nonlinear behavior around 400,000 ft 
starting altitude. Of note in this analysis is the very marginal 
instability or stability that may exist in this region. In all cases the 
real part of the eigenvalues are very near the imaginary axis when in this 
altitude region. In the case of the constant thrust rocket starting at 
400,000 ft, the periodic branch was found to have a region of stability. 
This implies that trajectories within a relatively small domain of 
attraction would be drawn to the limit cycle therefore the vehicle could 
expect to experience stable periodic motion, on the order of the orbital 
period, with fairly significant amplitude for velocity and altitude. 
Looking at the behavior of the translational states and given that the 
period of the limit cycle is nearly equal to that of a circular orbit, it 
seems likely that the limit cycles associated with the "resonance 
altitude" describe the velocity and altitude of an elliptical orbit. 

The bifurcation sweeps using the throttle parameter showed for the 
most part that the throttle is a very benign way to control the energy of 
the vehicle. Little was found that was of physical interest from the 
stand point of examining nonlinear behavior. The most interesting 
behavior was found for the air-breathing case with starting altitudes 
below 400,000 ft. At these "lower" altitudes the higher atmospheric 
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density made the aerodynamic forces more effective and the equilibrium 
energy management is done by changing velocity and altitude without 
rotating the vehicle. Only one physically realistic Hopf point was found 
for these low altitude sweeps. The periodic motion looked similar to that 
obtained from the body flap sweeps in terms of the altitude and velocity, 
however the rotational states experience a relatively high frequency 
sub-oscillations relative to what was seen before. In further comparing 
the results to the bifurcation sweeps using the throttle parameter and the 
body flap parameter for the constant thrust rocket case the pitch angle 
and the variation of the lift coefficient play key roles in the dynamics 
of the system. It seems without the pitch angle providing the impetus for 
instability the Hopf bifurcation at high altitudes where aerodynamic pitch 
damping is lost is not seen, which is not altogether surprising. 

Augmenting the vehicle model to obtain system stability for cases 
where the height mode is unstable seems to be intractable without changing 
the model to allow for commanded attitude input and the ability to 
generate or obtain the measurements of states or some value associated 
with a state that will allow for the minimization of some error. 

Certainly the rich dynamics associated with nonlinear phenomena has 
been demonstrated by the resulting complex behavior present even in this 
simple example. This work stands in contrast to those in previous studies 
who claimed to have explored the nonlinear nature of the longitudinal 
dynamics of a powered lifting hypersonic vehicle by simply including 
second order terms in the Taylor series expansions for small perturbation 
analysis. While no great departures of physical significance were found 
in this study from that which was previously obtained by others, this work 
does display many of the major findings of previous works and adds insight 
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to the expected local behavior of trajectories. 

Perhaps most useful of all, is this work displays the tremendous 
utility and encompassing nature of bifurcation analysis as well as the 
ease with which it can be applied to problems of this sort. 

Recommendations 

This work really stands as a first step. It opens the door to a 
variety of areas for further investigation. Several of these are briefly 
mentioned below. 

1. Add the lateral equations of motion and study the dynamics of a 
powered lifting hypersonic vehicle flying a minor circle. This would 
provide significant coupling of the longitudinal and lateral dynamics and 
should make for some interesting behavior. 

2. Define the aerodynamic forces and moment coefficients in nonlinear 
terms as found in references such as Etkin’s text (5:199,393), 

3. Include the rotation of the Earth in the equations of motion. 

4. Increase the accuracy of the thrust laws to reflect more up-to-date 
propulsion concepts such as ramjets and scramjets. This would bring Mach 
number and additional altitude dependencies. 

5. Develop higher order control systems to stabilize the height mode 
without destabilizing the phugoid mode. This will require the addition of 
states to the model to allow for commanded attitudes and feedback to 
controls with direct influence over altitude and velocity. Recall for 
this study the simple feedback of available states to the body flap proved 
worthless for stabilizing the height mode and of little value in 
maintaining phugoid stability much beyond an original Hopf point. 
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Appendix A: FORTRAN Listing of the User Supplied Subroutines for AUTO 


CURRENT AS OF 23 Nov 1990 


SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) 


This subroutine evaluates stationary solutions from the equations of 
motion for a powered lifting aerospace vehicle flying along a great 
circle about a nonrotating spherical Earth. 

Input parameters : 


NDIM 

U 


PAR 
PAR(1) 
ICP 
I JAC 


Dimension of U and F. 

State Vector containing U. 

V/VO the velocity along the flight path divided 
by a constant (nondimensional) 

Angle of Attack [radians] 

Pitch Rate of the Body relative 
to the Earth [rad/sec] 

Pitch Angle [radians] 

Radius from Earth's Center divided by a 
constant (nondimensional) 

Array of parameters in the differential equations, 
df Body Flap deflection [degrees] 

PAR(ICP(1)) is the initial 'free' parameter. 

*1 if the Jacobians DFDU and DFDP are to be returned 
-0 if only F(U,PAR) is to be returned in this call. 


U(1)« 

V/VO 


by a 

U( 2 )«* 

alpha 

U(3)~ 

q 

U( 4 )- 

theta 

U(5)« 

r/RO 


i 


Values to be returned : 

F - F(u, par) the right hand side of the ODE. 

DFDU - The derivative (Jacobian) with respect to U. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,K0,Ky,mslug,Lift 

DIMENSION U(NDIM),PAR(20) 

DIMENSION F(NDIM),DFDU(NDIM,NDIM),DFDP(NDIM,20) 

COMMON /CNST/ alt,L,Ky,K0 f re,R0,VO,gs,IRSTST,ITEST 
COMMON /FUNVL/ rho,g,Cd,Cl,CM,TPM,DR0DU5,DGDU5,S,mslug 
COMMON /AERO/ CdO,Cda,CIO,Cla,CMO,CMa,CMdf,CMq 
COMMON /CFORC/ orbper,Drag,Lift,altw,W,Fc,Tx,Ty 

. Set flag for Subroutine Const to use current values 

. of the States rather than initial values. 

ITEST-0 

. Call Subroutine CONST to obtain the States, plus the 

. necessary constants and functional values 

CALL CONST(U) 

************************************************************************** 

************** System of 5 Nonlinear Equations of Motion ***************** 
************************************************************************** 

. dv/dt SCALED ie U(l)- V/VO (Note U(l) is nondimensional) 

. NOTE: TPM - T/m 
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F( 1) ■={ (TPM) *DCOS (U{2)) - (0.5D0*rho*S*Cd* (V0*U( 1)) **2 )/mslug 
& -g*(DSIN<U(4)}*DC0S(0(2))-DC0S(U{4))*DSIN(U(2))) ) / VO 

.. d(alpha)/dt 

F{2)=U(3)+( g/(U(l)*V0) )*(DCOS(U(4))*DCOS(U(2 ;)+DSIN(U(4)) 

& *DSIN(U(2)))-(TPM/(U(1)*V0))*DSIN(U(2)'- 

& (0.5D0*rho*S*Cl*U(1)*V0/mslug) 

.. dq/dt 

F( 3 )«=( ( 0.5D0*rho*S*L*CM* (U( 1) * VO) * * 2) / (mslug*Ky**2) ) 

& -( (1.5D0)*( g/(U(5)*R0) )*(KO)*DSIN(2.0D0*U(4)) ) 

.. d(theta)/dt 

F(4)«U(3)+( ( (U(1)*V0)/(U(5)*R0) )*< DCOS(U(4))*DC0S(U(2))+ 
&DSIN(U(4))*DSIN(U(2)) ) ) 

.. dr/dt 

F(5)=((U{1)*V0)/R0)*(DSIN(U(4))*DC0S(U(2))- 
5 DCOS(U(4))*DSIN(U(2))) 

IF(IJAC.EQ.O)RETURN 


RETURN 

END 

SUBROUTINE STPNT(NDIK,U,PAR) 


.. In this subroutine the steady state starting point must be defined. 

.. (Used when not restarting from a previously computed solution). 

.. The problem parameters (PAR) may be initialized here or else in INIT. 

NDIM - Dimension of the system of equations. 

U - Vector of dimension NDIM. 

Upon return U should contain a steady state solution 
corresponding to the values assigned to PAR. 

PAR - Array of parameters in the differential equations. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,KO,Ky,mslug 

DIMENSION U(NDIM),PAR(20) 

COMMON /CNST/ alt,L,Ky,KO,re,RO,VO,gs,IRSTST,ITEST 
COMMON /FUNVL/ rho,g,Cd,Cl,CM,TPM,DR0DU5,DGDU5,S,mslug 
COMMON /AERO/ CdO,Cda,CIO,Cla,CMO,CMa,CMdf,CMq 

.. Initialize the problem parameters. 

PAR(1)-0,0D0 

write(*,*) 'Enter initial par(2)«Kq and par(3)~W/S 
& par(4)» Trho, and par(S)-throttle' 
read(* , 5 )PAR( 2) 
read(*,5) PAR( 3) 
read(*,5 )PAR( 4) 
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read(*,5)PAR(5) 
5 format(dl5.6) 

U( 2 )«=0.0D0 
•U(4)=0.0D0 
ITEST=10 
CALL CONST(U) 

C 

RETURN 

END 

C 

SUBROUTINE INIT 


. In this subroutine the user should set those constants that require 

. values that differ from the default values assigned in DFINIT. 

. (See the main documentation for the default assignments). 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20) 

COMMON /BLCDE/ NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT 
COMMON /BLDLS/ DS,DSMIN,DSMAX,IADS 
COMMON /BLLIM/ NMX,NUZR,RLO,RLl,AO,Al 
COMMON /BLMAX/ NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC 

★ ***************************a**;*A******************ifrdrri:*i fc ****ifc******* 

C******************* KEAD AUTO PARAMETERS *************************** 
C******************************************************************** 

C 

OPEN(UNIT=27,FILE-'DS.DAT',STATUS*'OLD') 

REWIND (27) 

C. ITEMS IN COMMON BLBCN - BASIC CONSTANTS 

READ( 27 ,*) NDIM 
READ(27,*) IPS 
READ(27, * ) IRS 
READ(27,*) ILP 
READ(27,*) ICP(l) 

write(*,*) 'Which parameter to vary?(PAR(2)«Kq PAR(3)»W/S 
£. PAR(4)-Trho, PAR(5)-Throttle setting)' 


read(*,5) I 
5 format(Il) 

ICP(2)=I 

C. ITEMS IN COMMON BLCDE - DISCRETIZATION CONSTANTS 


READ(27, *) NTST 
READ(27, *) NCOL 
READ(27,*) IAD 
READ(27,*) ISP 
READ(27,*) ISW 
READ(27,*) IPLT 

C. ITEMS IN COMMON BLDLS - STEPSIZE ALONG SOLN BRANCHES 

READ(27,* ) DS 
READ(27,*) DSMIN 
READ(27,*) DSMAX 
READ(27,*) IADS 

C. ITEMS IN COMMON BLLIM - LIMITS 

READ(27,*) NMX 
READ(27,*) NUZR 
READ(27,*) RLO 
READ(27,*) RLl 
READ(27,*) AO 
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READ(27,*) Al 

C. ITEMS IN COMMON BLMAX - MAXIMA 

READ(27,*) NPR 
READ(27,*) MXBF 
READ(27,*) IID 
READ(27,*) ITMX 
READ(27,*) ITNW 
READ(27,*) NWTN 
READ(27,*) JAC 
CLOSE (27) 

C 

RETURN 

END 

C 

FUNCTION USZR(I,NUZR,PAR) 


. This subroutine can be used to obtain plotting and restart data 

. at certain values of free parameters. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,Ky,KO,msiug 

DIMENSION U(5),PAR(20) 

COMMON /CNST/ alt,L,Ky,KO,re,R0,VO,gs,IRSTST,ITEST 
COMMON /FNVL/ rho,g,Cd,Cl,CM,TPM,DRODU5,DGDU5,S,mslug 

. Initially, for the steady state analysis, set NUZR-0 in INIT. 

. Then the functions specified below will be ignored. 

. When computing the branch of periodic solutions, set NUZR-4 in INIT. 

. Output will then be written in unit 8 for the values 

. of PAR(*) specified below. 

. Note that PAR(ll) is normally reserved. It is used by AUTO to keep 

.. track of the period (See main documentation). 

GOTO(1,2,3,4)I 

1 USZR*=PAR( 11) - 10.0 

RETURN 

2 USZR«PAR(11) - 14.0 

RETURN 

3 USZR=PAR(11) - 20.0 

RETURN 

4 USZR-PAR(ll) - 30.0 

RETURN 

END 

SUBROUTINE CONST(U) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,K0,Ky,Lift,mdot,m,mslug 
C 

DIMENSION U(NDIM) 

COMMON /BLBCN/ NDIM,IPS,IRS,ILP,ICP(20),PAR(20) 
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COMMON /CNST/ alt,L,Ky,K0,re,R0,V0,gs,IRSTST,ITEST 
COMMON /FUNVL/ rho,g,Cd,Cl,CM,TPM,DRODU5,DGDU5,S,mslug 
COMMON /AERO/ CdO,Cda,CIO,Cla,CMO,CMa,CMdf,CMq 
COMMON /CFORC/ orbper,Drag,Lift,altw,W,Fc,Tx,Ty 
C 

C************** Aerodynamic and Geometric Constants *********************** 
£******★*******************************************£******■***************** 
Cd0=0.0133D0 
Cda«=0.4D0 
C10=0.05D0 
Cla=0.5D0 
CM0=0.0D0 
CMa=-0.0548D0 
CMdf=CMa*l.5D0 
CMq=-0.028D0 
C 


C. Characteristic Length of Vehicle - overall length L=50 ft 

C 

L=50.0D0 

C.. Weight, Mass and Area of Vehicle 


C 

Wsl=700.0D3 


S“Ws1/PAR(3) 
mslug=Wsl/32.174D0 
C 

C. Radius of gyration in pitch [ft] - Ky~2 « Iy/m 

C 

Ky=25,0D0 
C 

C. K0°(Ix-Iz)/Iy 

C 

K0=-0.94D0 
C 

C. Radius of the Earth [ft] (standard geoid) 

C 

re-2.0903264468D7 

C 

C...... Gravity at Earth's surface [ft/sec‘2] 

C 

gs=32.174D0 
C 


V0«1000.0D0 
R0=1000.0D0 

raddeg»2.0D0*3.14159265359D0/360.0D0 
C 

C ************************************************************************* 

C ************ read altitude and set initial U(5) ************************* 
c ************************************************************************* 
c 

XF(XTEST .GT. 5) THEN 

OPEN(UNIT“25,FILE®'ALT.25',STATUS='OLD') 

REWIND (25) 

READ(25,*) alt 
C 

C. SET INITIAL VALUE FOR THE RADIUS U(5). NOTE: SCALED DOWN BY RO 

C 

U(5)«(re+alt)/R0 
C 
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************************************************************************* 
************** Density and Gravity variation with altitude *************** 

rt************************************************************************ 


IF( (U(5)*R0 - re) .LT. 6.0D5 ) THEN 
.. Convert altitude to SI units {km} 

altsi=( < U(5)*R0 - re ) / 3.2808399D0 ) / 1000.0D0 
.. Calculate the density in SI units {Kg/nT3] 

.. Constants used in polynomial (U.S. Standard Atmosphere Supp. 1966) 

co=o.ioooooooood+oi 

C1=0.3393495800D-1 
C2=-0.3433553057D-2 
C3=0.5497466428D-3 
C4 —0.3228358326D-4 
C5=0.1106617734D-5 
C6=-0.2291755793D-7 
C7=0.2902146443D-9 
C8—0.2230070938D-11 
C9=0.1010575266D-13 
CIO—0.2482089627D-16 
Cl1=0.2548769715D-19 

eqnl=CO+(Cl*altsi)+(C2*altsi**2)+(C3*altsi**3)+{C4*altsi**4)+ 
&(C5*altsi**5)+(C6*altsi**6)+(C7*altsi**7)+(C8*altsi**8)+ 

&(C9*altsi**9) + (C10*altsi**10)+(Cll*altsi**ll) 

rhoSI=l.2250D0/eqnl**4 

.. Convert from Kg/nT3 to Lbra/Ft“3 
.. (0.062427961 (lbm/ft A 3)/(Kg/jn*3) ) 

rho=rhoSI*0.062427961D0 

.. Convert to [Slugs/ft'3] 

rho=rho/32.174D0 

ELSE 

.. Equation for rho if altitude is greater than 600,000 ft 
rho= 2.16871253724D-10 * 

& DEXP{-8.89837671693D-6 * (U(5)*R0 - re) ) 


END IF 

. Calculate the gravitational acceleration 

g- gs * ( re / (U(5)*R0) )**2 

. Calculate the Atmospheric Pressure 
p0=2116.22D0 

* pi — 5.850746831820396D-05 

* p2=9.792179784448163D-01 

* p3=9.875326461241002D-05 
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p4=-6.044333173347913D-06 
p5=3.408653276857509D-09 
p6=-8.934489792146698D-07 
alt=(U{5)*R0 - re) 

Pa=pO*DEXP(pl*alt**p2)+p3*DEXP(p4*alt)+p5*DEXP(p6*alt) 

****************************************£****★**&*********************★*** 

******* Thrust Equations and Aerodynamic Coefficients ******************** 
************************************************************************** 

Cd=CdO+Cda*U(2)**2 
Cl=ClO+Cla*U(2) 

. Exhaust Nozzle Area [ft"2] 

An=40.0D0 

. Exhaust Velocity (Vexh) * 500 ft/sec 

Vexh=500.0D0 

. Rocket with constant Thrust at reference altitude and velocity. 

.Note: W/S-30 [lb/ft*2] at sea level, therefore 

. U(l)initial - sqrt{ (g r) / (1 + (rho r Cl S)/2m)) is 

IF (ITEST .gt. 5) THEN 

U(l)*= DSQRT( (g*U(5)*R0)/( 1.0D0+(rho*S*U(5)*R0*C1)/ 

& (2.0D0*mslug) ) ) 

U(3) - -U(1)/(U(5)*R0) 

. Thrust Constant 

T0=( 0.5D0*rho*S*Cd*U(l)**2 )/< rho**PAR(’4)) 
write(*,*) 'TO-',TO 

. Mass flow rate in (slugs/sec] 

. Assume the Pressure is expanded to the starting altitude value 

Pe*»Pa 

mdot=( (0.5D0*rho*(U(l)**2)*Cd*S) - (Pe-Pa)*An ) / Vexh 
write(*,*) 'MOOT*',indot 

. SCALE DOWN U(l) and nondimensionalize 

U(1)-U(1)/V0 

. PRINT RESTART DATA FILE 

OPEN(27, FILE**'REF.DAT' .STATUS**'NEW') 

WRITE(27,10) alt,TO 

1C FORMAT(40x,'REFERENCE VALUES',/,6X,'ALT [FT]', 

& 14X,'TO [ft‘4/sec*2]',/,2(IX,E15.8,4X)) 

CLOSE(27) 

OPEN( 40,FILE-'THRUST.DAT' ,STATUS**'NEW' ) 

OPEN(41,FILE-'COEFF.DAT',STATUS-'NEW') 

WRITE(40,21) 

WRITE(41,22) 

IRSTST-10 
try=0.OdO 

ELSEIF (IRSTST .LT. 5) THEN 
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OPEN(27,FILE='REF.DAT',STATUS-'OLD*) 

REWIND{27) 

READ(27,20) xjnkl,T0 

OPEN(40,FILE®'THRUST.DAT',STATUS®'NEW') 

OPEN(41,FILE='COEFF.DAT',STATUS*'NEW') 

WRITE(40,21) 

WRITE(41,22) 

20 FORMAT(/,/,2(1X,E15.8,4X)) 

21 FORMAT(6X,'H',15X,'Tx',15x,'D',15x,'Ty',15x,'L',15x,'W', 

& 15x,'Fc',14x,'Orb Per') 

22 FORMAT(8X,'H',15x,'Cl',15x,'Cd',15x,'CM',15x,'V',15x, 

& 'alpha',15x,'Rad') 

IRSTST=10 
END IF 
C 

Thrust = (T0*rho**PAR(4)) * PAR(5) 
orbper® (6.2831853072D0 * U( 5)*R0)/DSQRT(g*U(5)*R0) 

Drag® 0.5D0*rho*S*Cd*(U(l)*V0)**2 
Lift= 0.5D0*rho*S*Cl*(U(l)*V0)**2 
altw= U(5)*R0 - re 
w = mslug*g 

Fc = mslug*( ( (U(1)*V0)**2 ) / (U(5)*R0) ) 

Tx=Thrust*DCOS(U(2)) 

Ty=Thrust*DSIN(U(2)) 

TPM=Thrust/raslug 

if (try .It. l.OdO) write(*,*) 'TPM=',TPM 
try«=5.0d0 

. PAR(1) - body flap deflection 

. qO - Pitch rate due to spherical Earth in body axis system 

. Note all angles are in radians EXCEPT the body flap deflection 

. ( df=PAR(1) ) which is in degrees, therefore PAR(l) is multiplied 

. by 2*pi/360 [rad/deg) 

qO = -(U(1) * VO)/(U(5) * R0) 

U0 - DSQRT( g*U(5)*R0/((rho*U(5)*R0*Cl/2.0D0*mslug)+1.0D0)) 
df - PAR(1) - PAR(2)*(U(3)-qO) 

CM «* CM0+CMa*U(2)+CMq*(U(3)-q0)+CMdf*raddeg*df 

RETURN 
END 

subroutine BCND 
return 
end 

subroutine ICND 
return 
end 

subroutine FOPT 
return 
end 
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Appendix B: Standard Atmospheric Approximations for Density and Pressure 


The value of density for the Standard Atmosphere (14) is calculated using 
a different equation over two altitude regions. The lirst altitude range 
is from 0 to 600,000 ft. The density-altitude approximation for this 
range was obtained directly from the work by Vihn and Dobrzelecki (15:25) 
and provides values of density accurate to within 5% of the Standard 
Atmosphere for altitudes ranging from 0 to 200 Km (0 to 656,000 ft) 
(15:25). This equation is an inverse polynomial relationship given by: 


_P£J_ 

[A 0 + \Z + ... + A 11 Z 11 ] * 


(41) 


where: 


P 


Psl 


Z 


density [kg m’ 3 ] 

sea level density [kg m’ 3 ] 

altitude above standard geoid (6371.315 km) 

(note this is the average radius of the Earth at the 
equator, which is different for reference 15) 
Coefficients (j=l-ll) [km' 3 ] 


Table 2 

Coefficients for Density Polynomial 


0 

0.1000000000 

E 01 

1 

0.3383495800 

E-01 

2 

-0.3433553057 

E-02 

3 

0.5497466428 

E-03 

4 

-0.3228258326 

E-04 

5 

0.1106607734 

E-05 

6 

-0.2291755793 

E-07 

7 

0.2902146443 

E-09 

8 

-0.2230070938 

E-10 

9 

0.1010575266 

E-13 

10 

-0.2482089627 

E-16 

11 

0.2548769715 

E-19 


(15:26) 




r 


; Note the density is computed in [kg m" 3 ] and is then converted to [slugs 

ft' 3 ] using 1.9403232735 E-03 [{kg m' 3 )/(slugs ft" 3 )]. 

For altitudes above 600,000 ft the following exponential relation is used: 

p = 2.16871253724E-10 exp (-8.898376716935-06 Z) (42) 


where: p = density [slugs ft" 3 ] 

Z = altitude [ft] 

Note this equation yields the value of p in [slugs ft" 3 ] , 

Figure 44 shows the calculated density relative to the Standard 
Atmospheric data from reference 14. 

For the pressure altitude relation an exponential relationship was 
developed using the nonlinear fitting routine on MATLAB from MathWorks. 


P = PO ExpiPl Z P2 ) + P3 Exp(P4 Z) + P5 Exp { P6 Z) (43) 


where: P 

P0 
PI 
P2 
P3 
P4 
P5 
P6 
Z 


atmospheric pressure [lb ft ] 

2116.22 [lb ft' 2 ] 

-5.850746831820396 E-05 
9.792179784448163 E-01 
9.875326461241002 E-05 [lb ft" 2 ] 
-6.044333173347913 E-06 
3.408653276857509 E-09 {lb ft' 2 ] 
-8.934489792146698 E-07 

equatorial altitude above the Standard geode [ft] 


Figure 45 shows the quality of the pressure fit to the Standard 
Atmospheric data in reference 14. 
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Pressure [psf ] c 0ensity {3(ug/cu ft] 
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